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SUMMARY 


A method is presented for computing the wall coordinates of a hypersonic 
nozzle with real-gas^- effects. Results of calculations at a Mach number of 17 
for stagnation temperatures and pressures up to 5 >000° R and 1,000 atmospheres 
are presented. A procedure for calculating both the inviscid contour and 
boundary- layer displacement thickness is presented along with a complete com- 
puter program written in FORTRAN (FORmula TRANslation) language. Calculations 
are presented for a Mach number 17 nozzle for nitrogen at various stagnation con- 
ditions to indicate the difference between the use of real-gas properties and the 
ideal gas with constant heat-capacity ratio. The effect of stagnation conditions 
on both the inviscid flow field and the growth of the displacement thickness has 
been investigated. Whereas the present results were obtained for nitrogen, the 
method of calculation presented herein could be applied to other gases with only 
slight modification. 


INTRODUCTION 


The calculation of the inviscid contour of a hypersonic nozzle for a real 
gas involves an application of the method of characteristics in which the real- 
gas variation of the thermodynamic properties of the gas are considered. The 
physical-wall contour of the nozzle is then determined by adding the calculated 
displacement thickness to the inviscid contour. 

In reference 1 a procedure is presented for calculating only the inviscid 
contour of axi symmetric nozzles with reacting gases by application of the method 
of characteristics incorporating a variable isentropic exponent, whereas, in the 
present work the available thermodynamic data for an isentropic expansion are 
used directly in the computation scheme. The method for calculating the contour 
of axisymmetric nozzles for high-temperature air presented in reference 2 is 
similar to that presented herein but was developed independently. The differences 
between the present work and that of reference 2 include the manner in which the 

^The term "real gas" as used herein relates to the effects associated with 
high densities and also the variation of heat capacity with temperature. 


real-gas thermodynamic properties are employed and the detailed calculation pro- 
cedures for the boundary- layer determination. 

A comparison between the calculated inviscid coordinates of a hypersonic 
nozzle based on the real-gas thermodynamic properties of air at moderately high 
stagnation temperature and pressure and the coordinates based on ideal-gas rela- 
tions with constant heat capacity indicates that real-gas effects strongly affect 
the calculated nozzle contour. (See ref. 3 .) 

The procedure for calculating the wall contour of a hypersonic nozzle pre- 
sented in this report was used to design the nozzle of a Mach number 17 hyper- 
sonic facility which is presently being constructed for the Langley Research 
Center. This facility is to operate with nitrogen at stagnation temperatures 
up to 4,000° F and stagnation pressures up to 1,000 atmospheres for running 
times on the order of minutes and has a test- section diameter of approximately 
17 inches. 

It is the main objective of this report to present a method for calculating 
the physical-wall coordinates of a hypersonic nozzle which operates under condi- 
tions where real-gas effects are significant. In addition to the method of cal- 
culation, a comparison is given between the contours determined by inclusion of 
real-gas effects and of ideal-gas considerations with constant heat capacities 
for various stagnation conditions. The effect of wall temperature, size of inte- 
gration step, and other calculational restraints are also discussed. 

The method of calculation will be presented in enough detail to permit the 
direct use of this approach for determining the physical-wall contour of a hyper- 
sonic nozzle for a wide range of conditions. The modifications to this method 
which would be required for dealing with gases other than nitrogen and for con- 
siderably different conditions are also listed. 


SYMBOLS 


a velocity of sound 

B,C,D points in axisymmetric flow field as indicated in sketch (a); used as 
subscripts to indicate conditions at these points 

Cf skin-friction coefficient 

d* throat diameter 

H total enthalpy 

h static enthalpy 

M Mach number 

N exponent in velocity-profile relation 
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N Re local Reynolds number based on momentum thickness 0 

p pressure 

Q enthalpy correction term 

R gas constant for nitrogen 

r distance from source point in radial flow field 

r cr distance from source point to sonic sphere in radial flow 

r 1 radial distance from nozzle axis 

s distance along Mach line measured in meridian plane from nozzle axis 

T temperature 

T-j^ base temperature, 491.688° R 

T r gas temperature at nozzle throat (see eq. (Bl6)) 

u velocity 

u^ limiting velocity 

u T velocity in x direction 

W limiting velicity ratio, W = — 

U Z 

x distance along nozzle axis with x = 0 at nozzle throat 

x distance from source point of radial flow parallel to nozzle axis 

x distance measured parallel to nozzle wall with x = 0 at nozzle throat 

y distance perpendicular to nozzle axis 

y distance measured perpendicular to nozzle wall with y = 0 at nozzle 

wall 

d* 

y* nozzle throat height, — 

7 ratio of specific heats 

5 boundary-layer thickness 
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5* boundary-layer displacement thickness , defined by equation (B2) 

© characteristic temperature of molecular vibration for nitrogen, 

6,005-9° R 

0 momentum thickness 

flow angle 

p0 f i 

0 integrated values of flow angle, 0 = / ’ d0^, as given in 

1 1 

equation (a6) 

0* initial value of 0 at nozzle throat 

\i Mach angle 

p density 

\|r stream function 

\|r nondimen sional stream function, see equation (A9) 

cjo viscosity-temperature exponent 

Subscripts : 

aw adiabatic -wall condition 

L edge of laminar sublayer 

t stagnation conditions 

w wall condition 

1 inviscid free -stream conditions 
Superscript : 

o condition of low pressure 

METHOD OF CALCULATION 

The calculation of a real-gas nozzle contour is based on a characteristic 
solution for determining the inviscid flow field and then adding a correction of 
bound ary -layer displacement thickness to the inviscid nozzle contour. In the 
method presented herein, the actual thermodynamic properties of the expanding gas 


k 


are used to determine the inviscid flow-field boundary. A detailed procedure 
for calculating the inviscid portion of a hypersonic nozzle is presented in 
appendix A. 

After the inviscid contour has been determined by the procedure of appen- 
dix A, a displacement thickness based on real-gas properties is calculated from 
the edge of the inviscid contour to a physical wall. This calculation is based 
on a turbulent boundary-layer analysis in which the real-gas flow properties are 
used in a stepwise integration of the axi symmetric form of the momentum equation. 
(See pp. 393 to 395 of ref. 4.) The heat transfer to the nozzle wall is accounted 
for in the skin-friction law, and the skin-friction coefficient is obtained by a 
method presented by Persh. (See ref. 50 The quantities 5*/b and e/5 are 
determined by numerical integration with the use of a real-gas variation of den- 
sity through the boundary layer. The value of 0 at each point along the nozzle 
is obtained from the momentum equation by using an iteration scheme, wherein the 
first approximation to the radial distance is taken to be the sum of the inviscid 
radial coordinate, at the point of calculation, and the value of 5*, determined 
from the previous point of calculation. Successive iterations are made to deter- 
mine 5* within a set accuracy. A detailed description of the method of calcu- 
lating the displacement thickness is presented in appendix B. 

It is the aim of appendices A and B to present the methods of calculation in 
enough detail to permit one to calculate the coordinates of a hypersonic nozzle 
directly for the case of a real gas. 


RESULTS AND DISCUSSION 


Inviscid Results 

Attention is focused in this section on a number of calculated results which 
indicate how the various conditions and parameters affect the inviscid nozzle 
contour. First, a comparison was made between the calculated inviscid contour of 
a Mach number 17 nozzle based on real-gas properties for nitrogen and the contour 
based on the ideal-gas properties with a constant heat-capacity ratio of 7/5 for 
the same stagnation conditions. For this comparison, the stagnation pressure was 
taken as 1,000 atmospheres and the stagnation temperature was 4,200° R. The flow 
properties for the ideal-gas case were taken from reference 6, and for the real- 
gas case the thermodynamic data of references 7 &nd 8 were used. This comparison 
is shown in figure 1, where the throat height for both cases is equal to 
0.05 inch. It is noted that the height of the inviscid nozzle exit for the real- 
gas calculation is approximately 9 percent larger than the exit height based on 
the ideal-gas calculation. On the other hand, the contour height in the throat 
region is smaller for the real-gas case than for the corresponding ideal-gas case. 
This comparison indicates that the inviscid contours for the two cases at the 
same conditions and Mach number are significantly different, and it appears to be 
important to include the real-gas effects in hypersonic-nozzle calculations for 
stagnation conditions in this range. 
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Radial distance from nozzle axis, y, inches 


4.8 



Distance from throat along nozzle axis, x, inches 


Figure 1.- Comparison of inviscid contour as calculated for a real gas and an ideal gas at M-^ = 17> 

P t = 1,000 atmospheres, and = ^,200° R. 
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Two additional inviscid calculations were carried out in order to examine the 
effect of the choice of stagnation conditions on the inviscid contour of a Mach 
number 17 nozzle based on real-gas considerations, (it should be noted that the 
inviscid contour, based on an ideal gas with constant heat capacities, is inde- 
pendent of the stagnation conditions.) The effect of stagnation temperature was 
examined first. Figure 2 shows the inviscid contours for stagnation temperatures 
of 4,200° R and 5 >000° R, both for a stagnation pressure of 1,000 atmospheres. 

It is seen that the higher stagnation temperature gives a larger height at the 
nozzle exit and a decrease in the nozzle contour in the throat region. The effect 
of stagnation pressure was studied from a comparison of contours calculated for 
stagnation pressures of 3^0 atmospheres and 1,000 atmospheres, both at a stagna- 
tion temperature of 4,200° R. It is noted from figure 3 that the higher pressure 
results in a smaller inviscid exit height and a somewhat larger contour height in 
the throat region. 


Boundary -Layer Results 

In addition to the degree to which a nozzle contour is dependent on the real- 
gas effects in the inviscid region, the boundary- layer displacement thickness can 
also be influenced by real-gas effects within the boundary layer itself. The 
calculated boundary-layer displacement thickness will, therefore, be influenced 
by (l) the real-gas flow properties and their gradients along the edge of the 
inviscid flow field as determined from the inviscid calculations, and (2) the 
real -gas properties within the boundary layer. 

The degree to which these two effects influence the growth of the displace- 
ment thickness for a Mach number 17 nozzle with stagnation conditions of 
1,000 atmospheres and 4,200° R with a throat diameter of 0.10 inch was examined 
by comparison of the growth of 5* for three cases: (l) the inviscid flow field 
is calculated with the use of real-gas properties and the calculations within the 
boundary layer include the real-gas properties; (2) the inviscid flow field is 
again calculated with the use of real-gas properties, but the calculations within 
the boundary layer are based on an ideal gas with a constant heat-capacity ratio 
of 7/5; and (3) both the inviscid flow and the boundary layer are based on an 
ideal gas with a constant heat-capacity ratio of 7/5- This comparison is pre- 
sented in figure 4. It is noted from this figure that the largest difference is 
between the purely real-gas calculation and the purely ideal-gas result. At an 
axial distance from the throat of 90 inches, the value of 6* for the purely 
real-gas case is approximately 13 percent larger than for the purely ideal-gas 
case. On the other hand, figure 4 shows that the value of 5* at the same loca- 
tion for a real-gas boundary layer is only about 4 percent larger than that calcu- 
lated for the ideal-gas boundary layer with a heat-capacity ratio of 7/5, when 
both results are based on the same real-gas inviscid flow field. Based on these 
limited calculations, it appears that the real-gas effects in the inviscid flow 
field have a stronger influence on the determination of the displacement thickness 
than the real-gas effects within the boundary layer itself. 

An attempt was made to see why the displacement thickness is only slightly 
different for the calculation based on real-gas properties within the boundary 
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Radial distance from nozzle axis, y, inches 



Figure 2.- Effect of 


stagnation temperature on inviscid nozzle contour as calculated for a real-gas 
nozzle at Mj_ = 17 and = 1,000 atmospheres. 
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Radial distance from nozzle axis, y, inches 


5.2 



40 50 60 70 80 

Distance from throat along nozzle axis, x, inches 


120 


Figure 3.- Effect of stagnation pressure on inviscid nozzle contour as calculated for a real-gas 

nozzle at = 17 and T-j. = ^,200° R. 
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Figure 4.- Growth of houndary-layer displacement thickness as calculated for: (l) a real-gas inviscid 

contour with a real-gas boundary layer; (2) a real-gas inviscid contour with an ideal-gas boundary 

layer; and (3) an ideal inviscid contour with an ideal-gas boundary layer. p t = 1,000 atmospheres; 

T t = 4,200° R; T w = 6^0° R; d* = 0.10 inch. 


layer compared with the case for which ideal-gas properties with constant heat- 
capacity ratio were used within the boundary layer. The quantity 


1 T + 

2 



1 


which represents the ratio of the local static enthalpy, excluding the vibrational 
component, to the enthalpy including the vibrational mode was calculated across 
the boundary layer at several positions along the nozzle. The expression for the 
enthalpy, including the vibrational mode, is presented in reference 9- This 
quantity is unity when the vibrational component of enthalpy is zero and decreases 
continuously from unity as the vibrational component becomes more significant. 
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Figure 5 shows this quantity plotted as a function of the nondimensional 
distance through the boundary layer y/S at longitudinal distances from the 
nozzle throat of approximately 0.1, 1.0, and 100 inches. These calculations are 
based on a Mach number 17 nozzle with stagnation conditions of 1,000 atmospheres 
and stagnation temperature of 4,200° R and with a wall temperature of 650° R. It 
can be seen from figure 5 that the vibrational component of the local static 
enthalpy is always less than 10 percent of the total local static enthalpy 


Z T + 

2 



at any position within the boundary layer, even very near the 


nozzle throat. At a distance of 1 inch downstream from the throat, the enthalpy 
in the vibrational mode is less than 4 percent of the local static enthalpy and 
further downstream becomes even less significant. Because there is only a small 
percentage of the local static enthalpy at each point within the boundary layer 
in the vibrational mode, the small effect noted between the real -gas and ideal- 
gas boundary-layer results should be expected. The effect of pressure on the 
local static enthalpy within the boundary layer for the conditions of this example 
were found to be very small so that essentially all of the real-gas effects 
within the boundary layer are due to the vibrational component of enthalpy. 



Distance normal to wall, y/6 
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0 .2 .4 .6 .8 1.0 

Distance normal to wall, y/6 



0 .2 .4 .6 .8 1.0 


Distance normal to wall, y/6 


(a) x = 0.1 inch. (b) x = 1.0 inch. (c) x = 100 inches. 

Figure 5 .- Effect of molecular vibration on enthalpy through boundary layer at three locations ^ 
of x as calculated for a real -gas nozzle at M-^ = 17* P-^ = 1*000 atmospheres, T^. = 4,200 R, 

T w = 650° R, and d* = 0.10 inch. 
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In figure 6, a comparison of the calculated displacement thickness similar 
to that presented in figure 4 is given , the only difference being that the stag- 
nation temperature for the results shown in figure 6 is 5 >000° R. The comparison 
shown in figure 6 for 5^000° R is qualitatively the same as that shown in figure 4 
for 4,200° R. The comparison presented in figure 6 indicates, however, that the 
higher stagnation temperature results in a greater difference between the dis- 
placement thickness based on real-gas properties in the boundary layer as com- 
pared to those based on ideal-gas properties. This greater difference is believed 
to be due to the larger percentage of the local static enthalpy in the vibrational 
mode within the boundary layer at higher temperatures. 

In order to obtain an indication of the effect of stagnation conditions on 
the growth of the displacement thickness in a Mach number 17 nozzle, two sets of 



Distance from throat along nozzle axis, x, inches 


Figure 6.- Growth of boundary-layer displacement thickness as calculated for: (l) a real-gas inviscid 

contour with a real-gas boundary layer; (2) a real-gas inviscid contour with an ideal-gas boundary 

layer; and (3) an ideal inviscid contour with an ideal-gas boundary layer. p t = 1,000 atmospheres; 

T t = 5,000° R; T w = 650° Rj and d* = 0.10 inch. 
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calculations were made. In the first set, the stagnation pressure was chosen to 
he 1,000 atmospheres, and the stagnation temperature was set equal to 4,200° R 
and 5,000° R. For the second set of calculations, the stagnation temperature was 
chosen to be 4,200° R and the stagnation pressure was set equal to 340 atmospheres 
and 1,000 atmospheres. The results of these two sets of calculations are pre- 
sented in figures 7 and 8. It can be seen that an increase in stagnation temper- 
ature or a decrease in stagnation pressure results in a larger displacement thick- 
ness at each position along the nozzle. It should be remembered that these calcu- 
lations of 6* are also influenced by the changes in the inviscid flow field due 
to changes in stagnation conditions. 



Figure 7*- Effect of stagnation temperature on growth of real-gas boundary-layer displacement thickness 
as calculated for a real-gas nozzle at M]_ = 17 , p t = 1,000 atmospheres, T w = 650° R, and 
d* = 0.10 inch. 
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Displacement thickness, 6*, inches 



0 10 20 30 40 50 60 70 80 90 100 110 120 


Distance from throat along nozzle axis, x, inches 


Figure 8.- Effect of stagnation pressure oti growth of re&l-gas boundary- layer displacement as 

calculated for a real-gas nozzle at M]_ = 17, = 4,200° R, = 6^0° R, and d* = 0.10 inch. 

Because the physical -wall contour is the ultimate objective in the calcula- 
tion of a nozzle, figure 9 was prepared to show the difference between the 
physical -wall contours based on real-gas properties and on properties based on 
ideal gas with constant heat-capacity ratio of 7/5- This comparison is for a 
Mach number 17 nozzle with stagnation conditions of 1,000 atmospheres and 
^,200° R. It can be seen that the contour based on the ideal gas with constant 
heat-capacity ratio is significantly different from that based on the real-gas 
properties. 

The calculated values of 5* along the nozzle are based on a method that 
requires a number of assumptions which lead to limited certainty as to the correct 
local displacement thickness. Inasmuch as the physical contour is directly 
dependent on the calculated values of 5*, this same absolute uncertainty exists 
at each location along the physical contour. It should be noted, however, that 
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Distance from throat along nozzle axis, x, inches 


Figure 9*- Comparison of wall contour as calculated for an ideal-gas nozzle (ideal inviscid and ideal 
boundary layer) and a real-gas nozzle (real inviscid and real boundary layer) at M]_ = 17 , 
p^. = 1,000 atmospheres, = 4,200° R, T w = 6^0° R, and d* = 0.10 inch. 


the absolute values of 5* are rather small in the throat region so that the 
physical contour in the throat region is essentially determined only by the 
inviscid calculations. Further downstream, however, the boundary layer becomes 
a significant fraction of the physical contour. 

Effect of nozzle wall temperature .- The results of calculations presented up 
to this point are based on a constant wall temperature of 650° R. The effect of 
wall temperature on the calculated values of 5* was examined by comparison of 
results based on constant wall temperatures of 6 50° R and 1,500° R for a Mach 
number 17 nozzle with stagnation conditions of 1,000 atmospheres and 4,200° R. 
This comparison is presented in figure 10. It can be seen that an increase in 
wall temperature from 650° R to 1,500° R causes only a slight decrease in the 
displacement thickness along the nozzle. At the nozzle exit the difference is 
less than 1.0 percent, whereas, in the throat region the displacement thickness 
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Distance from throat along nozzle axis, x, inches 

Figure 10.- Effect of wall temperature on growth of boundary -layer displacement thickness as calculated 
for a real-gas nozzle at = 17, p t = 1,000 atmospheres, T t = 4,200° R, and d* = 0.10 inch. 

for the wall temperature of 650 ° R was less than 0.0002 inch smaller than the 
case for 1,500° R. Based on this calculation, the wall temperature does not 
affect the displacement thickness to a large extent. 

The results show that the displacement thickness was only slightly affected 
by a change in wall temperature from 650 ° R to 1,500° R; the boundary-layer pro- 
files of mass flow per unit area P uI / p l u l an( ^ P/P 1 for these two temperatures 

at 0 . 1 , 1 . 0 , and 100 inches from the throat were calculated as well as values of 
5*/S, 0/5, and 0 at the same nozzle positions. The profiles are presented 

in figures 11 , 12 , and 13, and the corresponding values of 5*/8, 0/&, and 0 

are shown in table I. It is noted from these three figures that profiles for 
both wall temperatures change appreciably along the nozzle but are qualitatively 
the same for both wall temperatures at a given nozzle position. The tabulated 
results of table I show that 5*/0 increases with wall temperature, but there 
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1.8 


1.0 



(a) Density ratio. (b) Mass rate of flow per area. 

Figure 11.- Boundary- layer profiles of density ratio and ratio of mass flow per area as calculated 
for wall temperatures of 6^0° R and 1,500° R at x = 0.1 inch, p t = 1,000 atmospheres, 

T t = ^,200° R, N = 7.58, and d* = 0.10 inch. 
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0 .2 .4 .6 .8 1.0 

Distance normal to wall, y/6 



0 .2 .4 .6 .8 1.0 

Distance normal to wall, y/6 


(a) Density ratio. (b) Mass rate of flow per area. 

Figure 12.- Boundary- layer profiles of density ratio and ratio of mass flow per area as calculated 
for wall temperatures of 650° R and 1,500° R at x = 1.0 inch, p t = 1,000 atmospheres, 

T t = 4,200° R, N = 6. 93, and d* = 0.10 inch. 
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(a) Density ratio. 


(b) Mass rate of flow per area. 


Figure 13.- Boundary-layer profiles of density ratio and ratio of mass flow per area as calculated 
for wall temperatures of 650° R and 1,500° R at x = 100 inches, p^. = 1,000 atmospheres, 

T t = 4,200° R, N = 5.83, and d* = 0.10 inch. 
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TABLE I 


BOUNDARY-LAYER PARAMETERS FOR A MACH NUMBER 17 NOZZLE AT 
P t x = 1,000 ATMOSPHERES AND T t ± = 4,200° R 


x, in. 

( approx . ) 

T w , 

°R 

6* 

6 

0 

6 

9 

e, 

in. 

0.1 

650 

0.065 

0.103 

0.63 

0.00095 


1,500 

.085 

•099 

.86 

.00090 

1.0 

650 

0.24 

0.079 

3-00 

0.0020 


1,500 

.28 

.074 

3.74 

.0016 

100 

650 

0 . 8 l 

0.012 

66.80 

0.044 


1,500 

.83 

.011 

79-00 

.036 


is also a simultaneous decrease in 0 which is assumed to he caused by a corre- 
sponding change in the value of Cf/2 determined from the skin-friction relation. 

The two effects tend to make 5* insensitive to wall temperature. 

The results of the work presented in reference 2 indicate that the choice of 
the skin-friction relation which is used in the boundary- layer calculation has an 
effect on the resulting values of 5*. It is shown in reference 2 that a change 
in wall temperature from 583° R to 1,500° R resulted in an increase in 5* when 
one skin-friction relation was used, whereas another skin-friction relation indi- 
cated a decrease in S* for the same temperature change. 

The skin-friction law of Persh (ref. 5) was selected for this present work 
based on satisfactory calibration results presented in reference 10. The nozzles 
tested in reference 10 showed good agreement with the designed performance based 
on Persh T s skin-friction relation for the calculation of the displacement 
thickness. 

In order that the nozzle wall temperature be more realistic in the actual 
nozzle operating conditions, a variable wall temperature was assumed along the 
nozzle with a hot wall in the throat region and a cooler wall downstream. The 
three wall-temperature variations which were used are plotted in figure ik and 
are labeled (a), (b) , and (c). Each of these temperature distributions begins 
at the nozzle throat with a temperature of 2,820° R and decreases continuously 
according to the equation indicated on the figure to asymptotic values for 
curves (a), (b), and (c) of 2,350° R, 1,500° R, and 650° R, respectively, further 
downstream. These wall-temperature variations are discussed in more detail in 
appendix B. The effect which these three wall-temperature variations have on the 
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Wall temperature, 



Figure lA.- Assumed variations of wall temperature with diameter ratio. 
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calculated displacement thickness along a nozzle is shown in figure 15* It is 



0 10 20 30 40 50 60 70 80 90 100 110 120 


Distance from throat along nozzle axis, x, inches 

Figure 15.- Effect of variable wall temperatures on growth of boundary-layer displacement thickness as 
calculated for a real -gas nozzle at M]_ = 17 , Pt = 1>000 atmospheres, T-j. = 4,200° R, and 
d* = 0.10 inch. (See fig. l4 for wall temperatures.) 


noted that the displacement thickness is only slightly influenced by the wall 
temperature. The reason for realizing only a small effect is believed to be due 
to the opposing changes in 6*/e and 9 with a change in wall temperature as 
discussed before. 

Effect of pressure on enthalpy .- As mentioned before, the real-gas effects 
considered in this work include the deviations from ideal-gas behavior due to 
excitation of the vibrational energy modes and to the high-density or high- 
pressure effects. Both of these effects are accounted for in the calculation of 
the inviscid flow field by use of the real-gas properties of nitrogen. In the 
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boundary-layer calculations , however, the vibrational effect on enthalpy is 
accounted for by an analytical expression with a function only of temperature. 

This analytical expression is then multiplied by a correction term which is the 
ratio of the actual enthalpy, including pressure effects, divided by the enthalpy 
based on the temperature alone. This pressure-correction term for the enthalpy 
is denoted as Q and is of the order of unity. The application of this correc- 
tion term is discussed in appendix B. The influence which the real-gas effects 
due to high pressure alone have on the calculation of 6* is shown in figure l6 
for a particular nozzle in the region near the throat. It can be seen that the 
difference between the result when using Q equal to unity (that is, no real-gas 
pressure effect considered) and the result when Q was taken into account, is at 
most a few percent for the conditions of this calculation. For higher stagnation 
pressures at the same temperature or for lower stagnation temperatures at the same 
pressure, that is, higher densities, the influence of using the actual values 
of Q would be greater. 


Computer Program 

The calculated results were obtained by use of an IBM 7090 electronic data 
processing system at the Langley Research Center. A listing of the detailed pro- 
gram statements in FORTRAN language (ref. 11) are presented in appendix C with 
appropriate comments for the benefit of agencies having access to digital com- 
puting machines. This program was set up and used to calculate hypersonic nozzle 
contours with nitrogen for stagnation temperatures up to 5 , 000° R and stagnation 
pressures up to 1,000 atmospheres. It should be pointed out that the method for 
calculating the inviscid flow field is general and can also be applied directly 
to other gases for stagnation conditions such that equilibrium dissociation 
effects must be considered. This method can, of course, also be applied to the 
calculation of the inviscid region of a hypersonic nozzle using helium at high 
stagnation pressures, for which case real-gas effects due to high density can be 
quite significant. The inputs to this present program for such cases would be 
the same as used in this work for the inviscid region, namely, the real-gas rela- 
tion between the Mach number and limiting-velocity ratio and the relation between 
the ratio of free- stream to stagnation density, and the limiting -velocity ratio 
for an isentropic expansion from a given stagnation condition. 

That part of the program which is concerned with the calculation of the 
boundary-layer displacement thickness, however, would in general require some 
modifications for other gases and much different stagnation conditions considered 
herein. A relation between the local static enthalpy within the boundary layer 
and the local pressure and temperature would be required which would replace equa- 
tions (B9) of appendix B. Accordingly, this would change the form of equa- 
tion (B12). The skin-friction law given by equations (B13), (Bl4) , and (B15) , 
could also be replaced by another law, if desired, or modified for a particular 
case. 
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Displacement thickness, 6*, inches 



0 .2 .4 .6 .8 1.0 1.2 1.4 1.6 1.8 2.0 


Distance from throat along nozzle sods, x, inches 

Figure 1 6 .- Effect of Q on growth of boundary-layer displacement thickness as calculated at M-j_ = 17, 
p t = 1,000 atmospheres, T t = 4,200° R, T w = 6^0° R, and d* = 0.10 inch. 
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CONCLUDING REMARKS 


A method for calculating a hypersonic-nozzle contour for a real gas has 
been developed and incorporated into a computer program which facilitates the 
rapid calculation of the inviscid flow field and the wall-boundary-layer displace- 
ment thickness. The procedure for calculating both the inviscid region and the 
displacement thickness is presented in enough detail to permit the direct appli- 
cation of this method. A working FORTRAN program for use on an IBM 7090 elec- 
tronic data processing system is also presented. 

The inviscid contour and the corresponding displacement thickness have been 
calculated for a number of conditions for nitrogen with stagnation temperatures 
up to 5,000° R and stagnation pressures up to 1,000 atmospheres for a Mach number 
of 17. The real-gas effects considered in these calculations are those that are 
associated with high-density gases and the variation of heat capacity with tem- 
perature due to vibrational excitation. This method of calculation presented 
herein may be applied to a number of systems with some modifications. Rather^ 
simple modifications are required for the consideration of equilibrium dissociated 
flow. 


Based on a number of calculations and the use of the present method, the 
following conclusions are indicated: 

1. A comparison between a calculated inviscid contour of a Mach number 17 
nozzle based on real-gas properties for nitrogen and on ideal-gas relations with 
a constant heat-capacity ratio of 7/5 for stagnation conditions of 1,000 atmos- 
pheres and 4,200° R shows that the inviscid contour based on real-gas considera- 
tions is considerably different from that found for the ideal gas. For this case 
the exit diameter of the inviscid nozzle is approximately 9-P ercen 't larger for 
the real-gas calculation than for the results based on ideal-gas relations. On 
the other hand, the diameter in the region near the throat is less for the 
inviscid real-gas result than for the ideal-gas result. 

2. It was found that when real-gas effects are taken into account the choice 
of stagnation conditions strongly affected the inviscid nozzle contour for the 
conditions examined. 

5 . The effect of stagnation conditions on the boundary-layer calculations, 
including real-gas effects, exhibited the same trends as for an ideal gas in that 
a decrease in the stagnation temperature or an increase in the stagnation pres- 
sure resulted in a decrease in the displacement thickness along the nozzle 
contour . 

4. The influence of the real-gas effects within the boundary layer on the 
calculation of the displacement thickness was much less significant than the 
influence of real-gas effects in the inviscid flow field on the determination of 
the inviscid nozzle contour. 

5. The calc\ilated displacement thickness for the condition of constant wall 
temperature was found to be only weakly affected by the level of the wall 
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temperature. For example, the displacement thickness, calculated for a Mach num- 
ber 17 nozzle with stagnation conditions of 1,000 atmospheres and 4,200° R and 
throat diameter of 0.10 inch, based on a wall temperature of 650° R was less than 
1 percent greater than the value for a wall temperature of 1,500° R at the nozzle 
exit. In the throat region the displacement thickness for the case in which the 
wall temperature was 650° R was less than 0.0002 inch smaller than the case for 
1,500° R. Also, calculations in which several realistic wall-temperature varia- 
tions along the nozzle contour were assumed indicated that the wall-temperature 
level and its variation along the nozzle do not strongly affect the calculated 
displacement thickness for the method presented herein. 


Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., January 7, 1965* 
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APPENDIX A 


CALCULATION OF INVISCID NOZZLE CONTOUR FOR A REAL GAS 


Characteristic Equations 

The procedure for calculating the inviscid nozzle contour for air as a real 
gas was briefly discussed in reference J. The approach used in this present work 
is essentially the same for the inviscid calculations. The basic equation for 
determining the inviscid contour for an axially symmetric nozzle is the potential- 
flow equation presented in reference 12 (p. 26l, eq. (^9))» The solution to this 
equation can be obtained through the method of characteristics which reduces to 
four characteristic equations which are readily adaptable to a finite difference 
technique. The derivation of the characteristic equations is also presented in 
reference 12. These four characteristic equations are: 


ay = 

dx 


tan(p +9) 


> first family 


^ - d0 tan n - h ^ = 0 
W y . 


(Al) 


and 


^ = tan(0 - p) 


„ second family 


— + d0 tan p 
W 



= 0 


(A2) 


where 


l 


1 " 


sin p sin 9 tan p 
cos(p + 0) 


and 


_ sin p sin 9 tan p 
2 cos(9 - p) 

These characteristic equations involve the five variables W, 0, x, y , and p , 
so that an additional relation is required. For the case of an ideal gas with a 
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constant specific-heat ratio, the expression relating 


M to W is 



and may he used because p is directly related to 


M by the equation 


(A3) 


p = sin -1 i- (A4) 

M 

The equivalent real-gas relation between W and M is obtained by taking the 
form of the ideal-gas expression as given in equation (A3) and tabulating the 
quantity l/w^ for various values of i/m^ as calculated for the isentropic 
expansion from a given stagnation condition. These are the values of l/w^ and 
l/M^ based on a real gas and are used directly in the calculation procedure with 
linear interpolation between the tabulated values. 


Determination of Flow Properties Along Flow-Region Boundaries 

The two families of characteristic equations (eqs. (Al) and (A2)) are solved 
by a stepwise procedure in which the initial values of the flow properties W, 

9, xj r cr , y/ T cr> M are first calculated along the edge of the character- 

istic flow fields. These initial flow properties are calculated along the line 
ABODE of sketch (a) . 


y/ r cr 



Sketch a 
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These quantities are based on real-gas considerations and are determined by a 
method presented in reference 13. For convenience in discussion, the nozzle is 
divided into four regions: I, first transition region j II, second transition 

region j III, radial flow; IV, uniform flow. Region III is bounded by the Mach 
lines BC and CD, and region IV is bounded by the final Mach line DE. The flow in 
region I is calculated by the method of characteristics and initial conditions 
along line ABC and the flow in region II is calculated by the method of character- 
istics and conditions along line CDE. The first task, then, is to calculate the 
flow conditions along line ABODE. 

Now the general equation for the radial-flow region as given by reference 13 
is 


(l _ M 2 )^u + 2 u = o 

V 'dr r 


(A5) 


For this region, the integrated characteristic equations, which are based on 
equation (A5 ) > are 


0 


I 



1 

2 


p W 

W(M=1) 


(m 2 - 1) 

w 


1/2 


dW 


(a6) 


and 



1 r w ( M 2 -_il aw 

2 ^w(M=l) w 


(A7) 


Because M and, therefore, W are known at point D (sketch (a)), the value of 
9 D at point D can be found by use of equation (A 6) and the real-gas relation 

between M and W. Values of W, 0|>, M, and r/r CI . along line CD are deter 

mined by choosing values of W which are successively less than the value of W 
at point D. Values of 0 f for particular values of W are calculated from equa- 
tion (A 6) and from the real-gas relation between M and W, along with the condi- 
tion that 9f = e i d " 9 I alon § tlie line CD. The corresponding values of r j r cr 

for these same values of W are determined from equation (A7) in a similar way. 
The coordinates for each of these calculations along line CD are determined from 
the equations 


-E— = — cos 0.p (A.8) 

r cr r cr 
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and 


= — — sin 9^ 
r cr r cr f 


(A9) 


The point indicated as C in sketch (a) is the nozzle -inflection point and 
corresponds to the maximum turning angle of the flow. The value of the flow 

angle at point C, 9 q is arbitrarily taken to be 12° in this work. It follows 

then that the integrated angle at point C , 9^- ^ , as obtained from equation (A 6), 

must equal the value determined from the condition that ®c = ®I D - ®I C’ 

values of M, W, xj r cr , and y/ v ct a ' t point C are fixed by equations (A7) , 

(a 8), and (A9). 

The flow properties along line DE are all constant and equal to conditions 
at point D and are determined by the procedure just mentioned. It follows that 
the quantities W, 9f, x/r cr , y/ r C r> M can 1:36 calculated along the 

boundary CDE and can serve as the starting conditions for calculation of flow in 
region II by the method of characteristics. 

The flow properties along line BC are determined by the same method as they 
were along line CD, except that the flow angle along line BC is obtained from the 
condition that 9-f = 9^- - 9-j- g, where p = ®I D - ®I C* 

The required flow conditions along line AB are established by assuming a 
linear Mach number distribution with respect to x/r cr from points A to B. This 

linear distribution is found by equating the slope of M with respect to x/r cr 
as determined at point B to the slope in this linear portion. The value of x/r cr 

for which the Mach number is unity is taken as the position of the throat of the 
nozzle. Inasmuch as the real-gas variation of W with respect to M has been 
determined, the establishment of a Mach number distribution along line AB also 
determines corresponding values of W along line AB. The necessary flow proper- 
ties along line ABC can, therefore, be calculated so that the method of character- 
istics can be applied to determine the flow field in region I. 


Characteristic Mesh Size 

The mesh size of the characteristic network is determined by the interval 
size between successive points chosen along the boundary line ABCDE. The inter- 
vals along lines BC and CD are determined by taking a constant interval of r jr cr . 

The interval size along line DE was taken to be the same as along line CD. The 
intervals between points along AB can also be chosen arbitrarily. For a given 
Mach number and set of stagnation conditions, the inviscid contour was calculated 
for various mesh sizes in which the interval in r / r cr along line ABCDE was set 

equal to 0.1, 0.2, 0.25, and 0-3. The effect of mesh size on the calculated 
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contour based on these four values was found to he negligible. In addition, the 
inviscid contour was calculated along line AB for intervals of r/r cr of 0.05> 

0.025, and 0.01 and for intervals along line BODE of 0.1. The inviscid contour 
was again essentially unaffected by this change in mesh size. It may be concluded 
that these mesh sizes are small enough not to influence the inviscid result. 


Calculation of Characteristic Network for Real Gas 

The characteristic network for a real gas is calculated by application of a 
method of successive approximations to equations (Al) and (A2) , along with the 
tabulated real-gas values of l/w 2 and l/M 2 and the flow conditions along line 
ABCDE. This method is outlined in pages 264-265 of reference 12. It has been 
determined by trial and error that three successive approximations are sufficient 
to give satisfactory results at each point of calculation within the character- 
istic network. 


Calculation of Streamline Along Inviscid Boundary 

Inasmuch as point C of sketch (a) is on the streamline which defines the 
inviscid boundary, the value of the stream function at any point along the invis- 
cid boundary must be equal to the stream function at point C. The procedure for 
determining the streamline which corresponds to the inviscid boundary, then, first 
involves the calculation of the stream function at point C. 

The differential form of the stream function in the radial-flow region as 
given in reference 15 is 


d\(f = pur 2 sin 0p d0f 


(A10) 


It is convenient to define a nondimensional stream function as 


f = 


t 


p t u Z r cr‘ 


(All) 


Inasmuch as the flow properties are constant for a given value of r/r QT in the 
radial-flow region, the differential form of the stream function may be integrated 
for a particular value of rj r cr to yield 


* = 



cos 



(A12) 


The values of W, r j T CT > and e f at P oirrt c can be found by a method already 
given. It is now required to know the real-gas relationship between pj p^_ 
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and W. For the case of an ideal gas with a constant heat-capacity ratio, the 
relation between P/P-^ and W is 

1 

£ - i 1 - w2 ) 7 ’ 1 <*») 

"G 

The corresponding real-gas relation between P/P^ anc ^ W is obtained from the 

form of equation (A13) and from tabulation of the quantity log p/p-^ for various 

values Qf log (l - W^) as calculated for the real-gas isentropic expansion from 
a given stagnation condition. A plot of log p/p-^ against log(l - W^) based 

on real-gas calculations showed that these quantities were nearly linear over 
small regions. Tabulated real-gas values of p/p^ and W were used with linear 

interpolation between log p/p^ and log(l - W^) for intermediate values of 
P/P^ and W for the calculation of \[ r at point C by equation (A12) . The value 

of the stream function \|r at any point on_the inviscid contour must equal the 
stream function at point C, that is, \|r = \|r^. 

The general differential form of the stream function for axisymmetric flow 
as given in reference 13 is 


d\|r = Pya ds 


(Al4) 


where y is the distance from the axis, a is the local speed of sound, and ds 

is the differential distance along a Mach line. This equation can be integrated 

along a Mach line where s is taken as 0 at the axis, so that the nondimen sional 

form of the stream function at any point s on a Mach line is 


♦ = 



(^)w(si n n) 



(A15) 


The value of \|r is calculated at each point of the characteristic network 
by stepwise integration of equation (A15) in which the real-gas relation between 
pjp ^ and W is used. The location of the intersection of the inviscid boundary 

with any Mach line is determined when \|r for a particular Mach line is equal 

to \|r . Points along the inviscid boundary are determined by interpolation 

L» 

between adjacent points in the network which have values of \| r that bound • 
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Real -Gas Thermodynamic Properties 


The calculations carried out in this work were based on the real-gas thermo- 
dynamic properties of nitrogen presented in references 7 and 8 . These data were 
plotted and the various thermodynamic quantities along an isentropic expansion 
were determined and tabulated. These tabulated quantities M, W, and p/p^ 

were then used directly in the machine calculations. The real-gas speed of sound 
for nitrogen above 100 atmospheres was based on the coefficients presented in 
reference 8. 
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APPENDIX B 


BOUNDARY-LAYER CALCULATIONS 


General Procedure 

After the inviscid contour has been determined from characteristic calcula- 
tions, the boundary- layer displacement thickness along the nozzle must be added 
to obtain a physical -wall contour. The coordinate system for calculating the 
displacement thickness 6* is shown in sketch (b) . 


r' 



Sketch b 


The flow quantities u^, p^, and at the edge of the inviscid region, 

and the flow angle 9^. are known from the inviscid results. The corresponding 
values of the static temperature at the inviscid boundary T^ was determined 
from tabulated real-gas values of T-^ against with linear interpolation 

between l/T-^ and M-^ for intermediate values. 


The momentum equation used in this work as derived in reference 4 is 


Cf 

2 


de 

dx 


(s ♦ 

\ 9/ u* 


du- 


1 + — 


dx 


dp 

p x d£ 


1 , 1 dr 1 
r dx 


(Bl) 
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where r' is the radial distance to the nozzle wall. The derivatives of veloc- 
ity, density, and radial distance with respect to distance along the nozzle 
wall x are calculated from differences between flow properties at successive 
points along the inviscid boundary. The values of p-^ and u-^ were taken from 

the inviscid results at various locations. 


Calculation of Shape Parameter for Real Gas 


The value of the shape parameter 6*/9 
the following integral equations: 


is calculated at each location with 


and 



(B2) 


(B3) 


The relation between the velocity profile and the variable of integration 
is assumed to be 



y/s 

(b*0 


The exponent N in this relation is obtained at each location along the nozzle 
from the following correlation: 


N = 1.77 log % e - 0.38 - |22. (B5) 

%e 

where N Re is based on local free-stream properties and the local value of the 

momentum thickness 9. This correlation is similar to that presented in refer- 
ence l4 and is based on experimental data for Mach numbers up to 9 • This corre- 
lation was used to determine values of N for values of N^ e that are beyond 

the data for which this correlation was determined; however, the values of 5*/5 
and 9/5 have been found to be quite insensitive to N at high values of N. 

The evaluation of equations (B2) and (B3) also requires a relationship 
between the and p/Pj_ within the boundary layer. This is obtained by 

beginning with the Crocco relation between total enthalpy and local velocity in 
the boundary layer as presented in reference 15 • 

H = A + (b6) 
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where the total enthalpy at any point in the boundary layer is 


H = h + (B7) 

The constants A and B in equation (B 6) are evaluated at the wall and at the 
edge of the boundary layer , so that the total enthalpy varies through the boundary 
layer according to the expression 


H - ^ + ( K i - 


(B8) 


where h^ is the static enthalpy at the wall. 

The expression for the local static enthalpy, which includes real-gas effects 
due to vibration but no real-gas effects due to pressure, can be written in non- 
dimensional form as 


h° 

7 T 

+ 

f © \ 

RT b 

2 T b 

V 

v e 0 / T - lj 


(B9) 


In order that the real-gas effects due to pressure on enthalpy be taken into 
account, equation (B9) is multiplied by a correction term which is equal to the 
actual real-gas enthalpy found in real-gas thermodynamic tables divided by the 
enthalpy calculated from equation (B9)> "both of which are taken at the same tem- 
perature. This correction term is denoted as Q = h/h°, where h° represents 
the enthalpy at low pressures. It follows that Q is a function of temperature 
and pressure or any two thermodynamic functions. Since the entropy is constant 
in the inviscid flow region, the quantity Q along the inviscid boundary can be 
written as a function of pressure alone for a given stagnation condition. This 
quantity can be represented by a power series in p^. 


— 1 • 0 t Ap^_ 


+ B 


(pi f + c (pi)' 


(BIO) 


The wall temperature is known, therefore the correction term along the nozzle 
wall can also be expressed in a power series in pressure: 


= 1.0 + Dp 1 + E (P ].) 2 + F ( P l ) 5 ( B11 ) 

where the coefficients A , B, C, D, E, and F in these equations are evalu- 
ated empirically from the real-gas thermodynamic data. 

The variation of Q through the boundary layer was taken to be linear with 
respect to h° from the wall value to the free-stream value Plots 
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of Q against h° for various pressures were found to be approximately linear 
over a rather large range. 

The equation which was used to relate the local velocity to the local density 
in the boundary layer was obtained by combining equations (B7 ) , (b8) , and (B9) 
and incorporating the pressure effect on enthalpy by the use of the function Q. 
The local temperature T within the boundary layer is replaced by the quantity 

T 

— -r — — so that the density and velocity within the boundary layer are related 

0/01 

by the equation 


Q 


7/W\ 



L u' \ + H 1 u' 

RT b\ " Vj RT b V 


1 K') 2 /u' f 

2 RTj) 


(B12) 


The solution to the quantities 5*/5 and e/5 as defined by equations (B2) and 
(B3) is obtained by a quadrature method at any value of x along the nozzle by 
simultaneous solution with equations (b 4) and (B12) . 


Skin-Friction Relation 


The skin-friction relation presented in reference 5 was used in this work to 
calculate the local skin-friction coefficient. This skin-friction relation is 
represented by the following three equations: 


Cf 

2 


= (20N) 


1-N 

1+N 


2 



T- 

T- 


N-ax>-i 
, N+l 


(B13) 


and 


t l . 

. T W 

T w - W u l\ 

^aw " -^li 


T 1 

_T 1 ‘ 

Ti Wy/ 

■ % 1 

W 




l-kjo 



Ur 


Un 


20N 


W R e(5/9) 


1 

N+l^V 


N+l 


(Bl4) 


(B15) 
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Equations (BlA) and (B15) are combined and solved by iteration at each point 
along the nozzle to determine T L . The resulting value of T L is then used in 

equation (B13) to calculate Cf/ 2 at that position. 


Calculation of Momentum Thickness 

The solution of equation (Bl) requires an initial value of the momentum 
thickness 0 at the point where x = 0. The effect of the magnitude of the 
assumed initial value of 0 on the growth of 0 along the nozzle was investiga- 
ted by choice of initial values of 0 of 10"^, 10"^, 10~\ and 10“^ inches. The 
results are shown in figure 17 which are based on the same inviscid conditions. 



Distance from throat along nozzle axis, x, inches 


Figure 17.- Growth of real-gas boundary- layer momentum thickness as calculated for several values of 
throat momentum thickness 0*, at M x = 17, p t = 1,000 atmospheres, T t = 4,200° R, T w = 6 50° R, 
and d* = 0.10 inch. 
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It is seen that initial values of 9 which are 10“^ or less do not significantly 
affect the values of 9 downstream from the throat. 

After the quantity d9/dx is determined from equation (Bl) at a particular 
location, values of 9 for adjacent downstream locations are found in a stepwise 
manner. 


Nozzle Wall Temperature 

The wall temperature enters into the calculation of 5* in equations (Bll) , 
(B12), and (Bl4). Calculations for a constant wall temperature are straight- 
forward. For a more realistic condition, the wall temperature was assumed to 
vary along the nozzle with the wall hot near the throat and cooler downstream. 

The expression used to describe these variations as a function of nozzle diameter 
ratio is 


T r - A(y /y - *) 1 - 8 
1 + B (y/y*) 1 * 8 


(Bl6) 


where T r is the gas temperature at the throat. The constants A and B are 

determined for an assumed wall temperature at the throat and an assumed wall tem- 
perature for a large value of y/y*, respectively. This relation given by equa- 
tion (Bl6) is based on a very rough heat balance between the nitrogen side of the 

nozzle and water-cooled passages in the nozzle wall. The term (y/y*) 1,8 results 
from the assumption that the heat-transfer coefficient on the nitrogen side of 
the nozzle is proportional to the area ratio raised to 0.9 power, as indicated in 
reference 1 6 . It should be noted that this expression for the variation of wall 
temperature along the nozzle was not developed to give an accurate prediction of 
wall temperature but, rather, to give a reasonable variation along the nozzle 
wall. 


Iterative Procedure for Determining Displacement Thickness 

The first estimate of 5* at a given nozzle position is found by assuming 
that the nozzle wall is represented by the inviscid contour. This initial value 
of 5* is then added to the inviscid contour to give a better approximation of 
the physical wall. The value of 5* based on this revised contour is then deter- 
mined and compared with the previous value. The various gradients are also 
revised in the course of this procedure. This iteration process is repeated until 
successive calculations of 5* at a given nozzle position differ by less than 
0.001 inch. 
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APPENDIX C 


COMPUTER PROGRAM FOR THE CALCULATION OF A NOZZLE CONTOUR 


A computer program was developed for calculating hypersonic nozzle contours. 
This program is written in FORTRAN language (ref. 11) and is presented at the end 
of this appendix. This program facilitates the computation of the inviscid noz- 
zle contour and the houndary-layer displacement thickness for a real gas based on 
the methods presented in appendixes A and B. 

Part I of the program determines the flow properties along the boundaries of 
the inviscid flow region. A tabulation of T, 1 /M^, l/w^, p/p^, based on a 

real gas, is supplied to the computer together with the flow angle at point C and 
the Mach number at point D (sketch (a)). The program prints the values of the 
flow properties of points A, B, C, and D. The program also writes on tape the 
flow properties along the boundaries DE, DC, BC, and BA. These properties are 
determined by the method described in appendix A. 

Part II of the program computes the inviscid nozzle contour. The tabulation 
of T, l/M^, l/w^, and p/p-f. supplied in part I is also used in part II. The 

tape written by part I and containing the flow properties along the boundaries is 
used by part II. Beginning at point D the method of characteristics described in 
appendix A is used to determine the flow properties along upward- sloping charac- 
teristic lines in region II. These lines are computed from the boundary DC and 
are extended -until the value of the stream function is equal to the stream func- 
tion at point C. The final points on these characteristic lines define the 
inviscid nozzle contour from point E to point C. Continuing from point B, char- 
acteristic lines are computed in region I from the BA boundary and are extended 
until the value of the stream function is equal to the stream function at point C. 
The final points on these characteristic lines define the inviscid nozzle contour 
from point C to the throat. The flow properties of the points defining the 
inviscid nozzle contour are printed and also punched on cards for use in part III. 

Part III computes the displacement thickness along the nozzle and applies it 
to the inviscid result to yield a physical contour. The flow properties at the 
edge of the inviscid region are supplied on cards punched by part II. The values 
of the constants used in equations (BIO) , (Bll) , (Bl6) and the variables such as 
p_k , u^, r^, cjd, T w , and T r are supplied to the computer. The values of the 

shape parameter S*/0, the skin-friction coefficient, and the momentum thick- 
ness 0 are calculated at each point on the boundary as described in appendix. B. 
An iterative procedure determines the value of the displacement thickness 8*. 

The value of 8* is added to the inviscid contour to give a better approximation 
of the physical wall. Finally, a second-order interpolation is employed to locate 
points on the physical wall at desired increments. These interpolated points on 
the wall are then printed. 

The following program has been used on the IBM 7090 electronic data process- 
ing system at the Langley Research Center to obtain the results presented herein. 
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THESE ARE INPUTS REQUIRED FOR PART I P-538.1 
E*TAB= 1/P2, INCREASING VALUES, LIMIT OF 1 OC VALUES 
WTAB*1/W2, INCREASING VALUES LIMIT OF 100 VALUES 
RHORT=RHO/RHOT, INCREASING VALUES LIMIT OF ICO VALUES 
EM A=MACH NO. AT POINT A T YP I C AL LY= 1 . 

THETC = THE T A AT POINT C IN RADIANS TYP I CALL Y= . 2094 395 1 
EMD=MACH NO. AT POINT D TYP ICALLY* 1 5. 

N*NUMBER OF VALUES IN TABLE OF 1/*2 

RL I M*L I M I T R/RCR ABOVE WHICH DELTA R/RCR IS KEPT 

WITHIN RERR OF DRORCR TYP I CALLY=4 .90 

RERR*ERROR LIMIT DR/RCR TYP I CALLY= .02 

DELW = DELT A W INTEGRATION INTERVAL TYPICALLY*. 00 1 

DRORCR =DEL T A R/RCR LIMIT 

ALSO DELTA R/RCR USED ON DE TYPICALLY=.2 
DDW = DE L T A W REDUCTION TYP IC AL L Y = .COOOO 1 
MAX = L I MI T TOTAL PCINTS=2000 

I BUG NOT = o FOR EXTRA PRINT, *0 OMIT EXTRA PRINT 

TTAB = T .DEGREES R, INCREASING VALUES, LI M IT OF 100 VALUES 

DXBA=DCLT A X ALONG PA TYP ICALLY=.G25 

DXDE*DEL T A X ALONG DE TYPICALLY*. 2 

TOTAL NO. POINTS ON DE* ( S/RCR* 10 ) /DXCE 


THESE ARE INPUTS REQUIRED FOR p AR T II P-538.2 
ARE A*G FOR CDE AREA, 1 . FOR ABC AREA, 2. TO END CASE 
START WITH CDE AREA 

PR I NT = 0 TO PRINT WALL POINTS 0NLY, = 1. TO PRINT NET AND WALL POINTS 
NAPRX=3=APPR0XI M ATI0NS PER POINT 
NN=NUMBER OF VALUES IN TABLE OF 1/M2 


TYPICALLY=6.3873E-7 
TYP IC ALLY* 1 . 

TYP ICALLY*180. 
TYPICALLY*. 00G1 


THESE ARE INPUTS REQUIREO FOR PART HI P-700.1 
RHOT *RHO SUB T TYPICALLY=6. 30045 1 
VE=V SUB L TYPICALLY=8977.4414 
R 1 =R SUB 1 TYPICALLY=55. 159 
Cl =WHERE MU=(C 1*T**3/2)/( T**D 1 +E 1 ) 

D 1 =CONST ANT IN MU EQUATION 
E 1 =CONST ANT IN MU EQUATION 
THI*THETA SUB I 
GC=32. 1739 
R= 1 774 .688 

THB*CHARACTERISTIC TEMPERATURE OF MOLECULAR VIBRATION TYP ICALLY =6005. 

OMEGA=EXPONENT IN V I SCOS l TY-TE “PERATURE RELATIONSHIP TYP I C ALLY* . 76 

T W = TW TYPICALLY = 650. R 

ERR=ERRCR TEST USE .00005 

XM=SCALE FACTOR TYP ICALLY*. 23832 

TL TG= 1 . =TL/T 1 TYP ICALLY* 1 . 

JL I M=TOT AL NO. OF POINTS ON WALL*NC. OF BINARY IN P UT CARDS 
DE BUG = 0 

N = 32N0 . OF POINTS PER INTERVAL GAUSS INTEGRATION 
L = 5N0. OF INTERVALS GAUSS INTEGRATION 
ACASE=CASE NUMBER 

ALPHA=WHERE Q=ALPHA+BETA«P+GAM*P**2+TAU*P**3 ICQ. BIO) TYPICALLY*!. 

BE T A = CONST ANT IN Q EQUATION T YP I CALLY* . 3546 3 1 25E-7 
GAM=CONS T ANT IN Q EQUATION TYP I C ALL Y=- . 1 53995 1 4E- 1 3 

T AU = CONST ANT IN Q EQUATION TYP ICALLY*. 42600556E-20 

PR=P SUB R TYPICALLY*. 72 

TT 1 G*T SUB T , 1 T YP I CALL Y = 54C 0 . 

TO*T SUB G TYPICALLY=491 .688 

AP=WHERE QW=AP+BP*P*CP*P*»2+DP*P**3 (EQ. BIO) TYP ICALLY* 1 . 

BP=CONST ANT IN QW EQUATION TYP ICALLY*. 18729E-6 

CP =CONST ANT IN QW EQUATION T Y P I CALL Y= . 1 1 277E- 1 2 

DP=CONSTANT IN QW EQUATION TYP I C ALLY*- . 1 4 72E-20 
HT=(H SUB T , 1 ) / I R*T SUB C) TYP I CALL Y =46 . 1 8 1 7 

APP=WHERE TW=(CPP-APPIY/Y*)**1.8)/(1+BPP(Y/Y«)**1.8) (EQ. B16) 

TYP ICALLY*0 

BPP = CONSTANT IN TW EQUATION TYP I CAL L Y = 0 

CPP = CONST ANT IN TW EQUATION=TR TYP I CALLY = 650 . 

YSTAR=CONSTANT IN TW EQUATION TYPICALLY* 1 . 

0PT = 0 IF C=Q I , NOT = 0 IF Q NOT*QI 
DIFN=ALLOWABLE DIFFERENCE IN N TYP I C ALL Y= . 0 1 
D I FC = ALLOW ABL E DIFFERENCE IN DELTA STAR RELATIVE TO 


TYPICALLY*. 0C1 


P-538. 1 

PART I COMPUTES POINTS ALONG DE , DC , BC , BA AND 
SAVES THEM ON TAPE9 

SUBROUTINE FINP IS FLOATING INPUT SUBROUTINE 
SUBROUTINE BIRD READS BINARY CARDS 
SUBROUTINE BIPUN PUNCHES BINARY CARDS 
SUBROUTINE GAUSS PERFORMS INTEGRATION 
SUBROUTINE RWCF FINDS SLOPE OF CURVE 

DIMENSION WTABl IOC ) , EMTABt 100) ,WLGGT( 100) ,RHOLNT( 100) , 


1 RHORT ( 100 ) , ALC ( 30) , ALW ( 40 ) . ALD ( 2 0 ) , W!2D0C) .EMC 2000) , 

2R0RCR! 2000) , THE T I ( 2000 ) , RORT ( 2000 ) , AK ( 2C0C ) , THET BC ( 1 5 ) ) , XDC ( 1 5? ) , 
3XBA( 1000) ,WBA( 1000) ,YDC( 150), WBC( 1 50 ) , THETDC ( 800 ) , XDC ( GOO ) , 

4 YDC ( 800 ) ,PSIOC!800),WDC(809),RORTOC(80U),YDE( *40 3 ) , XDC (3400) 

5 , T T AB ( 100),EMSQT( IOC), TEMP! IOC) 

DIMENSION DERI 110) 

EQUIVALENCE (EMTAB.DUMMYl ) , ( WT AB, 0UMMY2 ) , ( WLOGT, DUMMY 3 ) , ( RHOLNT , 
1DUMMY4 J , I TTAB,0UMMY5) , ( EMSOT , DUMMY6 ) 

EQUIVALENCE! W,XDE) , ( EM , XDE (2001) ) , ( FM ( 1 4f 1 ) ,YCE) , 

1 ( RORCR , YDE ( 60 1 )), (THETI .YDEI2601 ) ) , ( THET I (80 1 ) ,XBA) , 

2( THET I ( 1801 ) ,WBA) , ( RORT , WB A ( 20 1 ) ) 

COMMON EMTAB.WT A B, WLOGT, RHOLNT, TTAB.EMSQT 
COMMON RHORT, EMA, THETC , EMD.N.DSDE , DMBA , DW, DRORCR , 

1 ODW , MAX 

C READS INPUT CARDS 

C FINP IS FLOATING INPUT SUBROUTINE 

CALLF I NP ( 17, EMTAB.WTAB, RHORT, EMA, THETC, CMD.N.RL I M f 
1RERR,0ELW, DRORCR, DOW, MAX, I BUG , TT AB, DXB A , CXDE ) 

C WRITES HEAO 

WRITE OUTPUT TAPE 6,305 

C COMPUTES 1/T,M2,W,L0GE( 1-W2)»L0GE( RHO/RHC T ) 

DO 4 1=1, N 

TTABI I ) = 1./TTAB( I ) 

EMSOT! I1 = 1./EMTAB( I ) 

WN = SQRTF ( 1 ./WT AB ( I ) ) 

WLOGT! I ) =LOGF ( l.-WN**2) 

4 RHOLNT! I ) =LOGF ( RHORT ( I ) ) 

C REVERSES 1/T TABLE 

K = N 

DO 1000 1=1, N 
' TEMP(K)*TTAB( I ) 

1000 K=K— 1 

DO 2000 1=1, N 
2000 TTAB! I )=TEMP! I ) 

C REVERSES M2 TABLE 
K = N 

DO 3000 1*1, N 
TEMP(K)=EMSQT(I) 

3000 K=K- 1 

DO 4000 1*1, N 
4000 EMSQT! I ) =TEMP ( I ) 

C COMPUTES 1/MA2 AND INTERPOLATES 1/WA2 

0VMA=1 ./EMA**2 

C FTLUP IS FLOATING TABLE LOOK UP AND INTERPOLATION SUBROUTINE 

CALL FTLUPIOVMA.OVWA, 1 , N, EMTAD , WT AB ) 

C COMPUTES 1/MD2 AND INTERPOLATES 1/WC2 
0VMD=1 . /EMD* *2 

CALL FTLUP(OVMO,OVWD, 1 ,N, EMTAB.WTAB) 

C COMPUTES WA 

W A = SQRTF ( 1 ,/OVWA) 

C COMPUTES WD 

WD=SQRTF( l./OVWD) 

C SET R/RCR= 1 

RORCR! 1 )=1. 

C SET SUM*0 

SUM*0 

C SCT SUM2*0 

SUM2-0 

C INTEGRATION BY TRAPEZOIDAL RULE FROM WA TO WD OF (M2-U/W AND SQT 

C (M2-1J/W 

ZTERM=(EMA**2— l.J/WA 
ZTERM2=(SQRTF(EMA*»2-1.) )/WA 
DW*DELW 

C COMPUTE WA POINT 

W( 1 )=WA 
IC0UNT*0 

C TURN OFF LIGHT 

I F ( SENSE LIGHT 1)8,8 

C COMPUTES W,M,R/RCR,THETAI ,RHO/RHOT,K, LIMIT 0F2000 POINTS 

8 DO 12 1*2, MAX 

9 W( I ) =W ( I — 1 ) +DW 

I F ( W ( I )-WD ) 7, 7, 1 1 3 
113 I F ( SENSE LIGHT 1 ) 115,116 
1 16 DW=WD-W( 1-1) 

SENSE LIGHT 1 
GO TO 9 

7 OVW= 1 . /W ( I ) • • 2 

CALL FTLUPIOVW.OVM, 1 , N , WTAB , EMT AB ) 

EM(I)=SQRTFll./OVM) 

TEMP=(EM( I )*EM(I)-1. ) 

TEMP 1 =TEMP/W ( I ) 

TEMP2=SQRTF(TEMP)/W( I ) 

TEMP4= ( ZTERM +-TEMP1 )*.5*DW*SUM 
C EQUATION A7 

RORCR ( I)=EXPF( ,5»TEMP4) 

I F ( SENSE LIGHT 2)611,609 
611 SENSE LIGHT 2 

608 I F ( SENSE LIGHT 3)610,609 
610 RERR=RERR4.01 

I F ( SENSE LIGHT 2)609,609 

609 IFIRORCR! I )-RLIM) 10,602,692 

602 IF! ABSFI RORCR ( I )-RORCRI I - 1 ) -DRORCR )-RER«> ) 1C, 10,603 

603 I F ( RORCR ( I l-RORCRl I- 1 ) -DRORCR ) 604 ' 

605 DW=OW— DOW 

SENSE LIGHT 2 
IF(IBUG)5000, 5001, 5000 

5000 WRITE0UTPUTTAPE6,392,W( I ),DW, RORCR ( I ) 

5001 G0T09 

604 DW=0W4-DDW 
SENSE LIGHT 3 

IF! IBUG)5002, 5003, 5002 

5002 WRITE0UTPUTTAPE6,302,W! I ) , DW , RORCR ( I ) 


5003 G0TC9 

10 TCMP3=(ZTERM2+TEMP2)*.5*DW*SUM2 
I F ( SENSE LIGHT 2)506,606 

606 I F ( SENSE LIGHT 3)607,607 

607 I F ( IBUGI5C04, 5005, 5004 

50J4 WRITE OUTPUTTAPE6, 302, W( I ) , DW , RORCR ( I ) 

C EQUATION A6 

5005 THFTI ( I ) = • 5*TEMP3 
Z T ERM = TEMP 1 
ZTERM2=TEMP2 
SUM=TEMP4 
SUM2=TEMP3 

WLOG-LOGFl l.-Wl I)**2) 

CALL FTLUPIWLOG, RHOL N , 1 , N , WLOGT , RHOL NT ) 

RORTl I ) =EXPF ( RHOLN ) 

I COUNT = I COUNT* 1 
IFIIC0UNT-2000) 12, 12, 1600 

12 AK{ I ) = (CVW-1 . ) / OVM 

C I COUNT =N0 OF THETAS AND RORCRS STARTING WITH THE T I ( 2 ) AND R0RCRI2) 

C COMPUTE THE T AID ANO THETAIC 

115 THETID=THETI ( ICCUNT+ 1 ) 

C COMPUTE THETAI ,W,R/RCR, M ,K AT POINT B 

THETIB=THETID-2.*THETC 
THETIC=THETID-THETC 

CALL FTLUP(THETIB,WB, 1, I COUNT ,THETI(2),W(2) ) 

CALL FTLUP(THETIB,RORCRB, 1 , I COUNT , THET I ( 2 ) . RORCR I 2 ) ) 

0VWB=1 ./WB**2 

CALL FTLUP(OVWB,OVMB, 1 ,N, WTAB, EMTAB) 

EMB=SQRTF ( 1 . /OVMB ) 

AK B= ( OVWR- 1 • )/OVMr 

C FIT CURVE TO 7 POINTS AND READ OK/DMB AT POINT B 

C STORE KS IN ALO BLOCK 

J = 2 

DO 100 1 = 2,1 COUNT 

IFIEHI I*3)-EMB) 100,100,102 

102 ALD( J)=EM( I ) 

ALDl J - 1 ) * A K C I 1 
J* J*2 

IFIJ-8) 101,103,101 

103 J* J*2 

101 IFIJ-16) 100,104,104 
100 CONTINUE 

104 ALD ( 8 ) =EMB 
ALD l 7 ) *AK8 

C SHARE SUBROUTINE RWCF IS USED FOR SECOND DEGREE POLYNOMIAL 

C LEAST SQUARES CURVE FITTING ROUTINE USING ORTHOGONAL POLYNOMONI ALS 

CALL CF2F1 ( 0 , ALC , 0 , ALW , ALO , 2 , 7 ) 

CALL CF2F2(EMB,ALC,2,DER, 1 ,2) 

DKDMB=DER(3) 

C COMPUTE W,R/RCR,M, RHO/RHOT ,K,XBAR,Y,PSI AT POINT C 

CALL FTLUPI THETIC.WC, 1 , I COUNT , THE T I( 2 ) , W C 2 ) ) 

CALL FTLUPITHETIC.RORCRC, 1, 1 COUNT, THET I (2), RORCR I 2) ) 

OVWC= 1 . / WC * *2 

CALL FTLUPIOVWC.OVMC, I.N.WTAB.EMTAB) 

EMC=SQRTF( 1 . /OVMC ) 

WLOGC=LOGF (l.-WC**2) 

CALL FTLUP(WLOGC,RHOLNC, 1 , N , WLOGT , RHOLNT ) 

RORTC=EXPF (RHOLNC ) 

AKC = ( OVWC— 1 • ) /OVMC 
C EQUATION A 9 

YC=RORCRC*S INF (THETC) 

C EQUATION A12 

PSIC=R0RTC*WC»R0RCRC»«2»( I .-COSF( THETC I ) 

C COMPUTE DM/DRORCR AT POINT B 

DMDR*2./( ( WB/EMB )**2* ( 1 .-EMB*«2 ) *RORCRB* ( . 5» 

1DKDMB-AKB/EMB) ) 

C COMPUTE R/RCR AT POINT A 

RORCRA=RORCRB-(EMB-EMA)/DMDR 
C COMPUTE X AT POINT C 

XC=RORCRC»COSF (THETC ) -RORCR A 

NABPTS=0 

NBCPTS=0 

NDCPTS=0 

K* 1 

J*1 

4001 L = ICOUNT ♦ 1 
DO 300 1=2, L 

4002 IFITHETI ( I J-THETIB) 17,16,16 
16 IF I THETI (D-THETIC ) 13,14,14 

C COMPUTE THETA, Y,W,X AT POINTS ALONG BC LIMIT OF 150 POINTS 

13 THETBCI J)=THETI (I)-THETIB 
C EQUATION A 9 

YBC ( J ) *RORCR ( I )»SINF(THETBC(J) ) 

WBC ( J ) =W ( I ) 

XBC ( J ) = RORCR ( I )»COSF(THETBC(J) ) -RORCR A 

NBCPTS=NBCPTS*1 

I F ( J- 1 50 1 I 300, 1300, 1600 

C IF NO. OF POINTS EXCEEDS STORAGE , PR I NT OUT HERE 

C AND STOP 

1600 WRITE OUTPUT TAPE 6,500, J,K, I , ICOUNT 
CALL EXIT 
1300 J*J*1 

GO TO 300 

C COMPUTE THETA, X,Y,PSI,W, RHO/RHOT AT POINTS ALONG CD LIMIT OF 800 

C POINTS 

14 THETOCCK)*THETIC-THETI< f) 

XDC ( K ) =RORCR ( I ) »CCSF ( THETDC ( K ) J-RORCRA 
C EQUATION A 9 

YDC 1 K ) =RORCR l I ) *SINF l THETDC ( K ) ) 

C EQUATION A 1 2 

PSIDC(K)=RORT( I )*W{ I ) *RORCR( I ) ••2* ( 1 .-COSF ( THETDC ( K ) ) ) 

WDC(K) =W( I ) 


NDCPTS=NDCPTS*1 
RORTDC(K)=RORTl I ) 

IF (K-300) 1400,1400,1600 
1400 K = K* 1 

GO TO 300 
17 NAEPTS=NA8PTS*1 

300 CONTINUE 
400 NAfi = 0 

DX = . 

XB=RORCRB-RORCRA 
SAVE=RORT( I COUNT* 1 ) 

SAVE2=R0RCR1 ICOUNT+1 ) 

C COMPUTE POINTS ALONG BA LIMIT OF 10C0 
00221= 1 , MAX 
XB A ( I )=XB-DX 
IFIXBAI I ) 1600, 601,601 

600 XBA ( I ) =0 

601 TEM =XBA( I )*DMDR+1 . 

OVM = 1 ./ITEM ) »«2 

CALLFTLUPI OVM , OVW , 1 , N , EMTAB , WTAB ) 

WB A ( I ) =SQRTF ( 1 ./OVW) 

I F C 1-1000) 1500, 15C0. 1600 
1500 NAB=NAB*1 

IFIXBAI I ) )22»23,22 

22 DX=DX*DXBA 

C COMPUTE S/RCR AT POINT D 

23 SORCR=SQRTF{ 12.»PSIC)/ISAVE * ( W0/EMD**2 )) ) 

C ENDE=NUMBER OF POINTS ON DE LIMIT OF 340C 

ENDE=I SORCR* 10. )/DXOE 

EMUD = ARTNCFl 1 . /EMD , SQRTF I 1 . /EMD ) * *2 ) 

NDE=ENDE 

DSSMU=ORORCR*SINFIEMUD) 

DSCMU=URORCR«COSF I EMUD ) 

FLT= 1 . 

DO 211 1*1, NDE 

C COMPUTE X , Y OF POINTS ALONG DE 
XDE 1 1 ) =XDC l NDCPT S ) ♦FLT •DSCMU 
YDEI I )=FLT«DSSMU 
IF( 1-3400)21 1 ,21 1 , 1600 
211 FLT=FLT* 1 . 

RORCRD=S AVE2 

C WRITES DM/DRORCR AT POINT B 

WRITE OUTPUT TAPE 6,3J1,DMDR 
XA*0 
Y A = u 

THET AA = 0 

XB=RORCRB-RORCRA 

YB = u 

YD*., 

C WRITES HEADING 

C WRITES THETA, X,Y,W,M,R/RCR FOR A,B,C,D POINTS 

WRITE OUTPUT TAPE 6,303 

WRITE OUTPUT TAPE 6 , 3 1 C , X A , Y A , EMA , WA , RORCR A , THET A A 
WRITE OUTPUT TAPE 6,303 

WRITE OUTPUT TAPE 6 , 3 1 1 , XB , YB , EMB , WB , RORCRB , THCT I B 
WRITE OUTPUT TAPE 6,303 

WRITE OUTPUT TAPE 6 , 3 1 2 , XC , YC , EMC , WC , RORCRC , THET I C 
WRITE OUTPUT TAPE 6,303 

WRITE OUTPUT TAPE 6 , 3 1 3 , XDC I NDCPT S ) , YD, EMD , WD , RORCRD, THET I D 

303 FORMAT I 1H 1 0X 1 HX 1 6X 1 HY 1 6X 1 HM 1 6X 1 HW 1 4X 5HR/RCR8X6H THE T A I ) 

310 FORMAT I 3H A* 1 E 1 6 . 8 , 5E 1 6 . 8 ) 

311 FORMAT I 3H B=1E16.8,5E16.8) 

312 F0RMATI3H C* 1 E 16 . 8 , 5E 1 6 . 8 ) 

313 FORMAT I 3H D= 1 E 16 . 8 , 5E 1 6. 8 ) 

306 N=4*NDE 

C WRITE TAPE 9 SAVING INFORMATION ABOUT ALL POINTS ON AB,eC,CD,OE 

REWIND 9 

WRITE TAPE 9,PSIC,WD,N, (XDEILl , YDE ( L ) , L = 1 , NDE ) 

K=NDCPTS 

DO 405 1=1 ,NDCPTS 

404 WRITE TAPE 9 , ( XDC ( K ) , YDC ( K ) , THETCC I K ) , WDC ( K ) , PS I DC ( K ) , 

1 RORTOC ( K ) ) 

405 K=K— 1 

Z E RO*C 

WRITE TAPE 9, IZERO, ZERO, ZERO, ZERO, ZERO, ZERO) 

N = 4 • l NBCPTS* 1 ) 

WRITE TAPE 9 f PSIC v N,(XBC(I) , YBC 1 1 1 , THETBC I I ) , WBC (I),I = 1,N8CPTS) 
WRITE TAPE 9,XC,YC,THETC,WC 
DO 403 1=1, NAP 

403 WRITE TAPE 9, XB A ( I ) , WBM I ) 

WRITE TAPE 9, ZERO, ZERO 
REWIND 9 

C TABLES ARE SAVED IN COMMON STORAGE FOR NEXT CHAIN 

CALL CHAIN (2,B3) 

302 F0RMAT16E19.8) 

304 FORMAT I4E19.8) 

301 FORMAT ( 6H DMDR=1E19.8) 

305 F0RMATI9H1P-538. 1/) 

500 FORMAT ( 4 I 12) 

END 


C P-538.2 

C PART II NOZZLE CALCULATION 

C TABLES ARE ALREADY IN COMMON STORAGE FROM PREVIOUS CHAIN 

C PART II COMPUTES POINTS ALONG WALL FROM E TO C 

C AND FROM C TO THROAT, 

C PSI ALONG W ALL=PS I AT POINT C 


C PUNCHES X,Y,M,W,RHO/RHOT, THETA, T OF WALL POINTS 

DIMENSION STOR ( 4000 ), TEMPI 10J.EMTAB! 1uG),WTAB< 100),WLOGT( IOj), 

1 RHOLNT ( 100) 

2.TTA8I 100) ,EMSQT( 100) 

EQUIVALENCE ( EMTAB .DUMMY 1 ) ,( WT AB , DUMMY2 ),( WLOCT , DUMMY 3 ),( RHOLNT , 
1DUMMY4) , (TTAB.DUMMY5) , I EMSQT , DUMMY6 ) 

COMMON EMT AB, WTAe. WLOCT, RHOLNT, TTAB, EMSQT 

COMMON XA, YA, TH A , W A , F A , PS I A , XB , Y B , THB , WB , XC , YC , THC , WC , FC , P S I C , A MC 
1 AMUC.SMUC, CMUC,TMUC,PPC,Y AV, THAV.WAV, AMAV, AMU AV.SMUAV.CMUAV, TM'JAV 
2STHAV.YBV, THBV.WBV.AMBV, AMUBV, SMUBV.CMUBV, TMUBV , S THC V , XL I M , YL I * , 
3AMLIM.WLIM.PPL IM.THL IM, TLIM.PL IM.AMUL.SML IM.CMLIM.TMLIM.THPM, 
4TANP, T HMM, SINM, COSM, TANM.XATY, AC L,BCM, TEMP, WD.P.N, I, J, AREA, PRINT, 
5NAPRX.NN, ALNW, ALNP 
6.SINP.COSP.STOR 
FL IMFIXC,XA,P)=P*IXC-XA)*XA 
REWIND 9 

C TAPE 9 FROM CHAIN 1 CONTAINS 

C REC 1 PS I C , WD, 4NDE , X AND Y OF ALL PTS D+l TO E 

C REC 2 X,Y, THETA, W , PS I , RHO/RHOT OF EACH PT D TOC-1 

C 

C REC ( 2 + NOCP TS ) 0,0, 0,0, 0,0 

C REC ( 3+NDCPTS ) PS I C , 4 ( NBCPTS+ 1 ) , X , Y , THETA. W CF ALL PTS B+l TO C- 1 

C REC ( 4+-NDCPTS ) X,Y, THETA, W OF PT C 

C REC ( 5+NDCPTS ) X,W OF EACH PT B TO A* 1 

C 

C RECI5+NDCPTS+NABPTS) 0,0 

28 CALL FINP14, AREA, PRINT, NAPRX.NN) 

IFIAREA-1.) 41,41,40 

40 CALL CHAIN! 1,B3) 

C READ IN A LINE TO STOR 

C RE AO IN PT ON B LINE TO C PT 

41 WRITE OUTPUT TAPE 6,102 
I F ( AREA) 18,19,18 

18 RE AC TAPE 9 , PL I M , N , ( STOR ( J ) , J = 5 ,N ) 

K = N*4 

I = N ♦ 1 

READ TAPE 9 , l STOR l J ) , J= I , K ) 

XL I M*S TOR l N* 1 J 
YL I M = STOR I N*2 ) 

CALL AMCMU(ST0R(N*4) ,SMLIM,TMLIM,AMUL,AMLIM,CMLIM) 

THLIM=ST0R(N*3) 

WLIM=ST0R(N*4) 

CALL PPTI WLIM.PPLIM) 

EMSQ=AMLIM*AMLIM 

C SUBROUTINE FTLUP INTERPOLATES IN TABLE 

CALL FTLUPIEMSQ.TLIM, 1 , NN , EMSQT I 1 ) . TTAB Cl)) 

TL IM= 1 ./TL IM 

CALL BIPUNIXLIM.TLIM) 

WRITE OUTPUT T APE6 , 1 00 , XL IM , YL I M , AML I M , WL I M , PPL I M , THL I M , TL I M 
READ TAPE 9, I STCR I 1 ) , STOR l 4 ) ) 

STOR l 2 ) *0 . 

STOR I 3 ) *0 . 

GO TO 22 

C READ IN LINE DE AND A POINT FROM DC 

19 READ TAPE 9 , PL I M , WD, N , I STOR I J ) , STOR l J ♦ 1 ) , J* 1 , N ,4 ) 

READ TAPE 9 , XC , YC , THC , WC , PS IC , PPC 

DO 21 J* 1 , N , 4 

C THETA*0 AND W=WD FOR ALL PTS ON DE 
STOR l J *2 ) =0 . 

21 STOR I J *3 ) =WD 

20 I F I ARE A ) 22,23,22 

22 READ TAPE 9, XC, WC 

C Y» 3,THETA=0,F=0,PSI*0 FOR ALL PTS ON BA 

YC a C • 

THC=0. 

FC=*0. 

psic*o. 

GO TO 24 

C READ IN NEXT POINT FROM DC AND COMPUTE NEXT LINE 

23 READ TAPE 9 , XC » YC , THC , WC , PS I C , PPC 
CALL AMCMU I WC , SMUC , TMUC , AMUC , AMC , CMUC ) 

FC=PPC»WC»YC*SMUC 

24 IFIWC) 25,28,25 

C AFTER C IS COMPUTED MOVE STOR TO B AND MOVE C TO A AND TO STOR 

C WHEN PSI REACHES LIMIT COMPUTE END PT 

25 IF! AREA) 10, 11,10 

10 N = 1 

GO TO 7 

11 N=5 

GO TO 7 

6 N=N+4 

7 XB=STOR!N ) 

YB=ST0RIN*11 

THB=ST0RIN*2) 

WB = STOR I N+ 3 ) 

XA=XC 

Y A= YC 

THA=THC 

WA = WC 

F A = FC 

PS I A*PSIC 

I F I AREA ) 12,13,12 

12 STOR l N ) =XC 
STOR I N+ 1 ) =»YC 
STOR l N+2 ) =THC 
STOR l N* 3 ) =WC 
GO TO 14 

13 STOR I N-4 ) =XC 
STOR I N-3 )=YC 
STOR l N-2 ) =THC 
STCRIN-1 )=WC 

14 CALL GENPT 


I F ( PR I NT ) 27,26,27 

27 WRITE OUTPUT TAPE 6 , 1 00 , XA , YA , THA , WA , f A , PS I A , XR , YB , TH&, WB , XC , YC , 
1 AMC,WC,PSIC,PPC,FC,THC,AMUC 

C COMPUTE ALONG AN UPWARO SLOPING CHARACTERISTIC LINE UNTIL PSI 
C LIMIT IS REACHEO 

26 IF l PSIC-PL IM I 6,8,8 
8 I F ( ARE A ) 15, 16,15 

15 STCR ! N *4 ) =XC 
STOR ! N*5 ) =YC 
STOR ! N*6 ) *THC 
STCR ( N*7 ) *WC 
GO TO 17 

16 STOR! N ) = XC 
STOR! N* 1 1 = YC 
STOR ! N*2 ) *THC 
STCR l N*3 ) = WC 

17 P*(PLIM-PSIA)/(PSIC-PSIA) 

XL IM*FLIMF (XC,XA,P ) 

YLIM=FLIMF(YC,YA,P) 

THLIM=FLIMF(THC,THA,P) 

WLIM=FLIMF(WC,WA,P) 

CALL AMCMUIWL IM,SMLIM,TMLIM,AMUL,AMLIM,CML IM) 

CALL PPT (WLIM.PPLIM) 

EMSQ=AMLIM*AMLIM 

CALL FTLUPIEMSQ.TLIM, l.NN.EMSQT! 1 ) , T TAB ( 1 ) ) 

TL IM* 1 «/TL IM 

C PUNCH X,Y,M,W,RHO/RHOT, THETA, T 

CALL BIPUN(XLIM.TLIM) 

101 FORMAT! 1H 3E16.8) 

C PRINTS X, Y,M,W,RHO/RHOT, THETA, T CF WALL POINTS 

WRITE OUTPUT T A0£6 , 1 00 , XL I M , YL I M , AML I M , WL I M , PPL I M , THL I M , TL I * 

100 FORMAT ( 1 H 7E 1 6 . e ) 

102 FORMAT (1H 7X1 HX 16X 1HY 1 6X 1HM 1 6X 1 HW 1 1 X8HRH0/RH0T9X5HTHE T A 1 3X 1HT) 

GO TO 20 

END 


SUBROUTINE PSI (WX,PPX) 

C SUBROUTINE PSI COMPUTES PSI GIVEN W AND RHO/RHOT 

DIMENSION STOR! 4000), TEMP! 10) ,EMTAB( 1 30 1 , WTAB l 1 00 ) , WLOGT l 100) , 

1 RHOLNT ( 100) 

2.TTAB! 100) ,EMSQTl 100) 

EQUIVALENCE l EMT AB , DUMMY 1 ) ,! WTAB , DUMMY2 ),! WLOGT, DUMMY3 ),! RHOLNT , 
1DUMMY4), !TTAB,DUMMY5), !EMSQT,DUMMY6) 

COMMON EMTAB, WTAB, WLOGT , RHOLNT »TTAB,EMSQT 

COMMON XA,YA,THA,WA, FA, PS I A , XB , YB , THB , WR , XC , YC , THC , WC , FC , P S I C , A “C , 
1 AMUC,SMUC,CMUC,TMUC,PPC,YAV, THAV.WAV, AMAV, AMUAV, SMUAV,CMUAV, TMUAV, 
2STHAV,YBV, THBV , WBV , AMBV , AMURV , SMUBV , C M UBV , TMU3V, STHBV, XL ! M, YLI M , 
3AMLIM.WLIM.PPL IM,THLIM,TLIM,PLIM, AMUL,SMLIM,CMLIM,TML IM.THPM, 
4TANP, T HM M, SINM,COSM,TANM,XATY, AC L,BCM, TEMP, WD,P,N, I, J, AREA, PRINT, 
5NAPRX,NN,ALNW, ALNP.STOR 
6 , S I NP , COSP 
CALL PPT!WX,PPX) 

FC=PPX»WX»YC*SMUC 

PSIC*PSIA*. 5 *!FA*FC)*SCRTF! I XC-XA ) • f XC-XA ) ♦ { YC-YA ) • I YC- YA ) ) 

RETURN 

END 


SUBROUTINE PPT!WX,PPX) 

C SUBROUTINE PPT COMPUTES RHO/RHOT GIVEN W 

DIMENSION STOR l 4000) .TEMPI 10), EMTAB! 1 00 ) , WT AB l 1 00 ) , WLOGT l 1 0 D ) , 

1 RHOLNT ( 100) 

2, T T AB ( 100) .EMSQTl 100) 

EQUIVALENCE l EMT AB .DUMMY 1 ) , t WT AB , DUMMY2 ) ,! WLOGT , DUMMY3 ),! RHOLNT , 
1DUMMY4 ) , (TTAB.DUMMY5) , IEMSQT , DUMMY 6 ) 

COMMON E”T AB, WTAB, WLOGT, RHOLNT, T TAB, EMSCT 

COMMON XA, YA, THA, WA.FA.PS I A, XB.YB, THB, WB.XC.YC, THC, WC.FC.PS I C, A MC, 
1AMUC,SMUC,CMUC,TMUC,PPC,YAV,THAV,WAV,AMAV, AMUAV, smuav.cmuav, TMUAV, 
2STHAV, YBV,THBV,WBV,AMBV,AMUBV, S m UDV,CMURV,TMUBV,STHPV,XLIM,YLIM, 
3AMLIM, WLIM.PPL IM.THL IM, TLIM,PLIM,AMUL,SMLIM,CMLIM,TMLIM,THPf', 
4TANP, T HMM, SINM, COSM, TANM.XATY, AC L.BCM, TEMP, WD,P,N, I, J, AREA, PRINT, 
5NAPRX.NN, ALNW.ALNP.STCR 
6, SI NP, COSP 
ALNW = LOGF( 1 ,-WX*WX ) 

CALL F TLUP!ALNW,ALNP,1,NN, WLOGT l 1 ), RHOLNT! 1 ) ) 

PPX=EXPF!ALNP) 

RETURN 

END 


SUBROUTINE THM ! THX , T ANX , COSX , AMUX , THMX , SI NX ) 

C SUBROUTINE THM COMPUTES COS, SIN, TAN OF l THETA-MU ) GIVEN THETA AND 

C MU 

DIMENSION STOR ( 4000 ), TEMP! 1 C), EMTAB! 1 00 ) , WTAB ! 1 00 ) , WLOGT ! 1 0 ? ) . 

1 RHOLNT ( 100 ) 

2 .TTAB! 100 ) ,EMSQT! 100 ) 

EQUIVALENCE l EMTAB .DUMMY 1 ) ,! WT AB , DUMMY 2 ),! WLOGT , DUMMY 3 ),! RHOLNT , 
1 DUMMY 4 ) , !TTAB,CUMMY 5 ) , IEMSQT, DUMMY 6 ) 

COMMON EMT AB,WT AO, WLOGT, RHOLNT, TTAB.EMSQT 

COMMON XA, YA,THA,WA,FA,PSIA,XB,YB,THB,WB,XC,YC,THC,WC,FC,PSIC,AMC, 
1 AMUC.SMUC, CMUC.TMUC.PPC, YAV.THAV.WAV, AMAV, AMUAV, SMUAV.CMUAV, TMUAV, 
2 STHAV.YBV, THBV, WBV, AMBV, AMUBV.SMUBV, CMUB V ,TMUBV, STHPV,XLIM,YLI M , 
3 AMLIM,WLIM,PPLIM,THL IM.TLIM.PL IM,AMUL,SMLIM,CMLIM,TMLIM,THPM, 
4 TANP.THMM, SINM, COSM, TANM , XAT Y , AC L , BC M , TEMP , WU, P , N , I , J , AREA , PR IN T , 
5 NAPRX.NN.ALNW, ALNP.STOR 
6 , S I NP , COSP 
T HMX* THX-AMUX 
COSX*COSFl THMX) 


SINX=SINF {THMX» 
T ANX=S INX/COSX 
RETURN 
END 


C 

C 


SUBROUTINE THP ( THX , T ANX » COSX , AMUX , THPX , S I NX ) 

SUBROUTINE THP COMPUTES COS, SIN, TAN OF ( THET A*mu ) GIVEN THETA AND 
HU 

DIMENSION STORI4GOS) .TEMPI UJ ,EMTABI 1 00 ) » WTAB ( 1 Ov ) , WLOGT 1 10G ) , 

1 RHOLNT I 100) 

2 » T T AR ( 100) ,E M SQT( IOC ) 

EQUIVALENCE I CMT AB , CUMMY 1 ) , I WT AB , DU M MY2 ) , ( WLOGT , DU M MY l ) , (RHOLNT , 

1 DU V M Y4 ) , ( TTAB.DUMMY5) , I EMSQT ,PUM“Y6 ) 

COMMON EM TAB,WTAB, WLOGT, RHOLNT, TTAB.EMSCT 

COMMON X A , YA , THA , W A , F A , PS I A , XB , YB , THB , WB , XC , YC , THC , WC , FC , P S I C , A MC , 
1 AMUC.SMUC, CMUC,TMUC,PPC,Y AV.THAV.WAV, AM A V, AMU A V,SMU A V.CMUAV, TMUAV, 
2STHAV , YBV , THP.V , WBV , AMBV , AMUBV , SMUBV, CMURV , TMUHV , STHP.V , XL I U ,YLIM , 

3A M LIM,ftLIM,PPLI M »THLI w ,TLIM f PLlM f AMUL»SMLlM f C M Ll M ,TMLIM f THPM f 

UTANP, THMM, SINM.COSM, T ANM , X AT Y , ACL , BCM , TEMP , WD , P , N , I , J , ARE A , PR I N T , 
5NAPRX ,NN, ALNW, ALNP.STOR 
6, S INP, COSP 
THPX=THX*AMUX 
COSX=COSF(THPX) 

SINX=SINF(THPX) 

TANX*S INX/COSX 

RETURN 

END 


SUBROUTINE AMCMU I WX , SMUX , TMUX , AMUX , AMX , CMUX ) 

C SUBROUTINE AMCMU COMPUTES MU AND M GIVEN W 

DIMENSION STOR I 4000) .TEMPI 10 » , EMTAB I 1 00 ) , WTAB I IOC ) , WLOGT { 1 OC ) , 
1RH0LNTI 100) 

2.TTABI 100) ,E M SQTI 100) 

EQUIVALENCE I CMTAB, DUMMY 1 ) , I W TAR , DUMM Y2 ) , I WLOGT , DUMMY 3 ) , (RHOLNT , 

1 DU M MY4 ) , I TTAB.DUMMY5 ) , ( EMSQT ,rUM M Y6) 

COMMON EMTAB, WTAB, WLOGT, RHOLNT, TTAB.EMSCT 

COMMON XA,YA,THA,WA,FA,PSIA,XR,YB,THB,WB,XC,YC.THC,WC,FC,PSIC,AMC, 
1 AMUC.SMUC.CMUC, TMUC.PPC.Y A V, THAV.WAV, AM AV.AMUAV, SMUAV.CMUAV, TMUAV, 
2STHAV, YBV.THBV.WBV, AMBV, AMUBV, SMUBV.CMUBV, TMUBV.STHPV.XLIM.YL I M , 
3AMLIM,WLIM,PPLIM,THLIM,TLIM,PLIM,AMUl,SMLIM,CMLIM,TMLlM,THpv, 

4TANP.THMM, SINM.COSM, TANM.XATY, ACL.BCM, TEMP,WU,P,N, I ,J, AREA, PRINT, 
5NAPRX.NN, ALNW, ALNP.STOR 
6 , S I NP , COSP 
0W=1./(WX*WX) 

CALL FTLUPIOW.OM, 1 , NN , WT AB I 1 ) , EMT AB I 1 ) ) 

S M UX=SCRTF(OM) 

A M X= 1 . /SMUX 
CMUX=SQRTF( l.-OM) 

TMUX=SMUX/CMUX 
A M UX*ATANF (TMUX) 

RETURN 

END 


C 

C 

C 


SUBROUTINE GENPT 

SUBROUTINE GENPT COMPUTES GENERAL POINT BY THREC DIMENSIONAL 

I RROTAT IONAL FLOW EQUATIONS MAKING AS MANY APPROXIMATIONS AS 
DESIRED 

AVRGFI THA,THC)*.5«(THA+THC) 

D I MENS I ON STORI4000) .TEMPI 10),CMTABI !OG),WTAB(100)»WLOGT(10S), 


1 RHOLNT I 1 OO ) 

2.TTABI 100 ) .EMSQT ( 100) 

EQUIVALENCE I E MT AB .DUMMY 1 ) , (WTAB,DU M MY2),I WLOGT , DUMMY 3 1 , I RHOLNT , 
1DUMMY4) , I TTAB.DUMMY5) , I EMSQT , DUMMY6 ) 

COMMON EMTAB,WTAB,WLOGT,RHOLNT,T TAR, EMSQT 

CO M MON XA, YA,THA,WA,FA,PSIA,XR,YB,THB,WB,XC,YC,THC,WC,FC,PSIC,AMC, 

1 AMUC , SMUC.CMUC , TMUC , PPC , YAV , TH AV , WAV , AMAV , AMUAV, SMUAV , CMU AV , TMU AV , 
2STHAV, YBV , THBV, WBV , AMBV , AMUBV , SMURV , CMURV , TMUBV , S THBV , XL I M , YL I« , 

3AMLIM,WLIM,PPLIM,THL IM.TLIM.PL IM.AMUL.SML I M.CMLIM.TML I M, THP*', 

4 T ANP, THMM, SINM.COSM, T ANM , X AT Y , ACL , BCM , TE MP , WU , P , N , I , J , ARE A , PR IN T , 
5NAPRX.NN, ALNW, ALNP.STOR 
6, S INP. COSP 
DO 51=1 , NAPRX 
I F I I - 1 ) 4,3,4 
3 THAV=THA 
THBV= THB 
WAV=WA 


WBV=W 3 

3 C CALL AMCMUIWAV, SMUAV, TMUAV, AMUAV, AMAV, CMUAV) 

CALL AMCMU I WBV, SMUBV, TMUBV, AMUBV, AMBV, CMUBV ) 

CALL THP I THAV, T ANP, COSP, AMUAV, THPM.S INP) 

CALL THHI THDV.TANM.COSM, AMUBV, THMM, S I NM) 
STHAV=SINF(THAV) 

STHBV=SINF(THBV) 

XATY=XA*TANP-YA 

XC=(XATY*YB-XB*TAN M )/I TANP-TANM ) 

YC=XC*TANP-XATY 
AC L=SMUAV* TMUAV* S THA V/CCSP 
PCM=SMUBV*TMUBV*STHBV/COSM 
YAV=AVRGF (YA.YC) 

Y 8 V=AVRGF(Y 8 ,YC) 

I F ( YAV+YBV ) 1 , 2,1 
1 TEMP=-THA*TMUAV+IXC~XA) *ACL/YAV 

THC=(-KA-WAV*TEMP + WB+WBV*(THB*TMUBV*(XC-XB)*BCM/Y"V) )/ 
1 (WAV*TMUAV*WBV*TMUBV) 

WC=WA*WAV* (TMUAV*THC+TEMP) 

GO TO 5 


2 THC=(-WA*WB)/(2.*(WAV*TMUAV*WRV*TMUBV) ) 
WC=WA*WAV*2.*TMUAV»THC 
CO TO 5 

4 THAV=AVRGF(THA,THC) 

THBV=AVRGF(THB,THC) 

WAV=AVRGF(WA,WC) 

WBV=AVRGF ( WB , WC ) 

GO TO 30 

5 CONTINUE 

CALL AMCMU(WC,SMUC,TMUC,AMUC,AMC,CMUC) 
CALL PSI(WC.PPC) 

RETURN 

END 


C P-700.1 i 

C PART III COMPUTES BOUNDARY LAYER USING 

C WALL POINTS PUNCHED BY PART II ( 

DIMENSION FS ( 2 ) ,SUM(2) ,ANS(2) , X I 2000 ) , Y ( 200C ) ,W(200G) ,RRT(2U00) , 

1 THFC 2000) ,TI (2000) ,SVY( 2000) .SVDELTC 2000) .DELS! 2000) ,YY (200 j ) ) 

COMMON OPT 

COMMONTI ,RRI f TO , THB , QSUM, CW , EN , ERR ,U I , HW, R , HT ,CON 1 , C0N2, 

ICON, T,H,Q,FR,FPR,H3,HWS,K,RH0I ,01 , X , Y , W , RRT , THF , SVY.SVDE LT, DELS , YY l 

2 ,RH0T,VE,R1,C1,D1,E1,THI,GC, 

3 OMEGA, TW, X« , TLTG, JL I M, DE BUG , N , L , AC ASE , ALPHA, 

*»BETA,GAM,TAU,PR,TT1G,X1 ,Y1 ,EM1 ,W1 ,RRT1 ,THF1 ,TI 1, 

5 TLTI ,TT1 1, THETA, DXB, DU, 

60R , DY , PI ,EMU,RETH,TEMP,FT,FTPR,H1 , DSDEL , : 

7TH0EL.TAW, AL P , CW , C 2 , C 3, C4 , FT T , F TTPR, H2 , CF2 , DSTH , 

8DTDX, YPD, OEL I 

9 , DP , CP , BP , AP , HI 

1 WRITE OUTPUT TAPE 6,25 i 

25 FORMAT ( 1 H 6X5HACASE13X4HRH0T14X2HVE1 5X2HR 1 1 5X2HC 1 1 5X2HD 1 
115X2HE 1/7X3HTHI 1 5X2HGC 1 5X 1 HR 1 5X3HTHB 

21 3X5H0MEGA 1 4 X2HT W 1 4X 3HERR/8X2HXM 1 4X4HTLTG 
312X5HALPHA13X4HBETA13X3HGAM14X3HTAU15X2HPR/ 

47X4HTT1G14X2HT015X2HAP15X2HBP15X2HCP 1 5X2HDP 1 5X2HHT ) 

103 FORMAT ( 1H 7E 17.8) 

100 FORMAT ( 1HO6E17.8/6E17.0) 

CALLFINP(39,RH0T,VE,R1,C1,D1,E 1 , TH I , GC , R , THE , OMEGA , 

1 T W, ERR, XM, TLTG, JL I M, DEBUG, N,L, AC ASE, ALPHA, BETA, GAM, 

2TAU,PR,TT 1G,TC,AP,BP,CP,DP,HT, APP , BP P , CPP , YST AR , 

30PT , D I FN , D I FD ) 

WRITE0UTPUTTAPE6, 1 03 , AC ASE , RHOT , VE , R 1 , C 1 , D 1 , E 1 , TH I ,GC , R , THB , OME GA , 

1 TW, ERR, XM, TLTG, ALPHA, BETA, GAM, TAU,PR,TT1G 

2,T(. ,AP,BP,CP,OP,HT I 

26 FORMAT! 1H 8X 1 HX 1 6X 1HY 1 4X5HTHET A 1 3X4HCF/2 j 

1 14X1HN14X6HTH/DEL/5X7HY+DELST 1 1X5HDELST 1 3X3HDEL 

214X4 HKETH13X4HTLTI 1 OX 9HDEL ST/DEL ) | 

J = JL I M 

C SUBROUTINE BIRD RE AOS BINARY CAf\CS CONTAINING X , Y ,M , W , RHO/RHOT, 

C THETA, T 

C LIMIT OF 2000 CARDS 

C RE AC IN WALL POINTS EXIT TO THROAT 

C MULTIPLY X AND Y RY XM AND STORE THROAT TOEXIT 

1000 CALL BIRD (XI, TI 1) 

X( J)=X1*XM 
YI J)*Y1»XM 
W( J)«W1 
RRT l J ) =RRT 1 
THF ( J ) =THF 1 
TIC J)*TI 1 
SVY( J)=Y( J) 

50 SVDELT ( J ) =0 
J = J-l 

IF ( J )8 1,81, 1000 
81 TLTI=TLTG 
TT 1 1=TT1G 
THET A= TH I 

WRITE0UTPUTTAPE6.26 
SENSE LIGHT 1 

2 008 J J = 1 , J L I M 
K = J 

RHOI*RRT(J) *RHOT 
UI*WCJ)»VE 
PI-RUTK JMRHOI 

EMU=(C1*(TI(J)**1.5) )/( (TI ( J)**D1 )+E 1 ) 

C EOUATION BIO 

IF(PI-20000. )2C3,203,204 

203 01=1. 

G0T04 

204 QI=PI*(PI*(PI*(TAU)*GAM)+BETA)+ ALPHA 

4 I F ( THB/TT 1 1-24. 1120, 120, 101 j 

101 FT=VE*VE-CI*(7.*TT1 1*R) 

FTPR=-QI*(7.*R) 

GOT0 102 j 

120 TEMP=EXPF (THB/TT11 ) { 

FT=VE*VE-CI»(7.*TT1 1 «R* ( 2 . *R»THB ) / ( TEMP- 1 . ) ) ! 

FTPR=-0I»(7.*R+( ( ( THB*THB*2.»R)/( TT1 1 • TT 1 1 ) )»TEMp) 
l/( (TEMP-1. )•• 2) ) 

102 Hi =-FT /FTPR 
TT11=TT11+H1 

IF(ABSF(H1/TT1 1 )-ERR)5,5,4 
C EQUATION B 1 6 

5 IF (CPP)40C,401 ,400 
401 TW = TT 1 1 

G0TC4.2 

400 TW= ( CPP-APP* ( SVY ( J ) /YST AR ) *• 1 . 8 I / ( 1 . +RPP» ( SVY ( J ) / YST AR ) •* 1 .8 1 
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C EQUATION 011 

402 IFIPI-20000. >205,205,206 

205 QW = 1 • 

GOT0207 

206 QW=PI*(PI*(PI*(CP)*CP)*BP)+AP 

207 IF (THB/TW-24. 1 104, 104, 105 

105 HWS=I7.*TW)/(2.*T0) 

GOTO 1 06 

104 HWS*(7.*TWI/(2.*T0 1+THB/ 1 ( EXPF ( THB/TW ) - 1 • > • TO > 

106 IFITHB/TI ( J1-24. ) 107, 107, IOC 

108 HI = (7.*TI(J) 1 / ( 2 . * TO ) 

GOTO 1 09 

107 HI = ( 7- *T I ( J) )/(2.*T0)*TH8/( (EXPF (THR/TIl J) 1-1. )*T0) 

109 QSUM= ( QI-QW 1/ (HI— HWS ) 

I F ( THB/TW-24 . 1110,110,111 

111 HW=( I 7 .*TW 1 / ( 2.*TU 1 )*CW 
GOTO 1 12 

1 10 HW=( (7.*TW)/(2.*T0 )*THB/( ( EXPF ( THB/TW J - I . I * TO J )*QW 

112 RET=(RHOI*UI 1 / C EMU* 1 2 . 1 

503 RETH=RET*THETA 

C EOUATICN B5 

EN=1.77*.43429448*L0GF(RETH)-.38-20Q./RETH 
I F ( SENSEL IGHT 1)511,502 
502 IF (ABSFIEN-ENP 1 -01 FN 1500, 500,501 
511 SENSEL IGHT 1 
501 RR I = . 125 

C SUBROUTINE GAUSS INTEGRATION 5 INTERVALS AND 32 POINTS PERINTERVAL 

CALL GAUSSIN.L.O. , 1. , 2 , S , F S , SUM , ANS , 0 1 
C DELTA STAR/DELTA 

DSDEL= 1 • -ANS ( 1 1 
C THETA/DELTA 

THDEL* 1 .-CSOEL-ANS ( 2 J 

TAW=( PR**. 333333331* ITT 1 1-TIC J1J+TII J1 

ALP=( 1.*0MEGA)/(EN*1.) 

CW=TH/TI(J1 
C2=(TW-TAW1/TI (JJ 

C3=(20.«EN*THDEL/RETH)*»( l./(EN*l.) 1 
C4 = ( T AW-T I ( J 1 J/TI l J) 

7 FTT«CW-C2*C3*ITLTI**ALP)-C4*C3*C3*(TLTI**<2.*ALP) J-TLTI 
FTTPR=-ALP*C2*C3*( TLTI •*( ALP-1 . 1 1 -2 . *ALP *C4 *C 3*C3 • 

1 (TLTI**!2.*ALP-1. 1 1-1. 

H2*-FTT/FTTPR 
C EQUATION B 1 4 

C TL/TI 

TL TI =TLT I ♦H2 

IF (ABSFIH2/TLTI J-ERR16,6,7 
C EQUATION B 1 3 

C CF/2 

6 CF2=( l./C20.*ENJ J*C3*C3*( CTLTI J**( ( 1 . *2 . *OMEGA-EN 1 
1 / ( EN* 1 . 1 1 1 

C DELTA STAR/THETA 

DSTH=DSDEL/THDEL 
500 DELS! J)=OSTH*THETA 
YPO-SVYl JJ+OELSI J! 

DEL=DELSI J J/OSDEL 

WRITE0UTPUTTAPE6,100,X( J» ,Y( JJ , THE T A , CF2 , EN , THOEL , YPO , CEL S ( J 1 ,D EL , 
1RETH, TLTI , OSDE L 
YYC J I =Y * J) 

Y { J J =YPD 

IFIJ-l 1504,505,504 

505 DELSP=DELS( J) 

G0T0509 

504 DXn=(X( J1-XI J-1 1 1/COSFITHFI J-l 1 1 
DU=VE*(W( J1-WC J-1 J J 
DR=RHOT*(RRT ( JJ-RRTI J-l 1 1 

DY=Y( J J -Y ( J-1 1 
C D THETA/D X 

DTDX=CF2P-THETAP*( I ( 2. +DSTHP 1 *DU J / ( U I P*DXB ) *DR/ ( RHO I P*CXB 1 
1+0Y/(YIJ-1 1 *DXB J J 
THETA=THETAP+(OTDX)*DXB 
I F { SENSEL IGHT 11508,506 

506 I F ( ABSFI (DELS! J1-DELSP1/SVYI JJ J -0 IFD 1 509 , 509 , 508 

508 DELSP=DELS( J) 

E N P = E N 
G0T0503 

509 YCJ+1 1 = Y I J ♦ 1 l+DELSIJJ 
CF2P=CF2 
THETAP=THETA 
DSTHP=DSTH 

UIP=UI 
RHOI P=RHO I 

510 SENSEL IGHT 1 
30 CONTINUE 

1001 X I NT=0 

WRITE OUTPUT TAPE 6,514 

C USE SECOND ORDER INTERPOLATION TO FINO GIVEN POINTS 

C DX=.10 IF Y LESS THAN 1. 

C DX=.25 IF Y GREATER THAN 1. AND LESS THAN5 . 

C DX=. 50 IF Y GREATER THAN 5. 

1005 CALL FTLUPIXINT,YINT,2, JLIM,X,YY) 

WRITE OUTPUT TAPE 6, 5 1 3 , X I NT , Y INT 
IF(YINT-YY( JLIM) 11006,1,1 

1006 IFIY1NT-1. J 1002,1009, 1009 

1002 XINT=XINT+.1 
GO TO 1005 

1009 IF (YINT-5. 1 1003, 1004, 1004 

1003 XI NT =X INT*. 25 
GO TO 1005 

1004 X I NT=X INT* . 5 
GO TO 1005 

513 FORMAT ( 1H 2F12.4) 

514 FORMAT ( 1 H 6X1HX11X1HYJ 


END 


SUBROUTINEFOFX(S,FS) 

DIMENSION F S C 2 ) ,SUM(2) ,ANS(2) ,X 12000) , Y ( 200C ) , WI 2 30C ) ,RRT(2 0.) , 
1THFI2JC0) ,TI (2000) ,SVY(2G00) ,SVCELT( 2000) , CELS (200i'l ,YY(20G. ) 
COMMON CPT 

COMMON T I ,RRI , T J,THR,3SUM,QW,EN,ERR,UI ,HW,R ,HT , CON 1 ,C0N2, 
1CON,T»H,Q,FR,FPR,H3,HWS,K,RHOI ,31 ,X, Y.W. RRT , THr, SVY.SVCELr. tels , yy 
CON 1 * ( (UI*UI )/ (R*TO) ) 

CC\2 = S • • ( 1 ./EN) 

CCN=C0N1*(HW*( 1.-CON2)/CONHHT*CON2/CON1-.5*CON2»CON2) 

C ITCRATfc FOR RHO/KHOI 

11 T*TI(K)/RKI 

C IF EXPONENTIAL GREATER THAN 24. OMIT TERM 

C EQUATION B9 

IF I THR/T-24. )113,113,114 
1 14 H=(7.*T)/(2.»T0) 

GO TO 1152 

1 13 H=( 7.*T)/( 2.»T0)4THP/ ( (EXPFl THR/TJ-1 . ) *TC ) 

1152 IF(0PT)115, 1150, 115 

1150 0 = 0 1 

GO TO 1151 

1 15 0= ( QSUM ) • ( H-HWS ) *QW 

C ir EXPONENTIAL GREATER THAN 24. OMIT TERM 

1151 I F ( TH3*RRI/TI <K)-24. )116,116,117 

117 FR*C* ( 1 7 . • T I (K) )/(2.»T0»RRI ) )-CON 
FPR=C* I ( -7 ,*T I C K ) )/( 2.*T0*RRl*RRI ) ) 

GO TO HE 

116 FR*Q* ( ( 7. *T I ( K ) )/(2.»T0«RRI ) ♦THE/ ( ( EXPF ( THB»RR I /T I ( K ) ) - 1 . ) • TO > > 
1-CCN 

FPR=Q* ( f -7. •TICK) )/( 2.*T0*RRI*RRI )-( THB«THR*EXPF ( THB*RRI / 

1TI (K) ) )/(TC*TI(K)*(EXPF(TH8*RRl/TI(K) )-l.)».2) ) 

118 H3=-FR/FPR 

C EQUATION B 12 

RRl\’=RRl+H3 
IF(RRIN) 125, 125,126 

125 RR I =RR 1 • . 5 
GO TO 11 

126 RR I = RK IN 

IF(ABSF(H3/RRI J-ERR) 10,10,11 
10 FS( 1 )=RRI*(S**( l./EN) ) 

40 FS(2J = RRI*(S«*(2./EN) ) 

HOLC=S 
H0LC2 = F S ( 1 ) 

44 H0LD3*FS(2) 

IOC FORMAT ( 1 H 7E17.8) 

RETURN 

END 


REFERENCES 


1. Guentert , Eleanor Costilow, and Neumann, Harvey E. : Design of Axisymmetric 

Exhaust Nozzles by Method of Characteristics Incorporating a Variable Isen- 
tropic Exponent. NASA TR R-33 , 1959* 

2. Enkenhus, K. R., and Maher, E. F. : The Aerodynamic Design of Axisymmetric 

Nozzles for High-Temperature Air. NAVWEPS Rep. 7395; U.S. Naval Ord. Lab. 
(White Oak, Md.), Feb. 5; 1962. 

3. Erickson, Wayne D. , and Creekmore, Helen S. : A Study of Equilibrium Real-Gas 

Effects in Hypersonic Air Nozzles, Including Charts of Thermodynamic Proper- 
ties for Equilibrium Air. NASA TN D-231, 19^0. 

4. Young, A. D. : Boundary Layers. Vol. I of Modern Developments in Fluid 

Dynamics - High Speed Flow, L. Howarth, ed. , The Clarendon Press (Oxford), 

1953; PP- 375-^75- 

5. Persh, Jerome, and Lee, Roland: A Method for Calculating Turbulent Boundary 

Layer Development in Supersonic and Hypersonic Nozzles Including the Effects 
of Heat Transfer. NAVORD Rep. 4200 (Aeroballistic Res. Rep. 320) , U.S. 

Naval Ord. Lab. (White Oak, Md.), June 7, 1956. 

6 . Ames Research Staff: Equations, Tables, and Charts for Compressible Flow. 

NACA Rep. 1135; 1953- (Supersedes NACA TN 1428.) 

7. Hilsenrath, Joseph, Beckett, Charles W. , et al. : Tables of Thermal Properties 

of Gases. NBS Cir. 564, U.S. Dept. Commerce, 1955; PP- 297-368. 

8 . Hall, N. A., and Ibele, W. E. : Thermodynamic Properties of Air, Nitrogen and 

Oxygen as Imperfect Gases. Tech. Paper No. 85 , Eng. Exp. Station, Univ. of 
Minnesota, Dec. 1951* 

9 . Donaldson, Coleman duP. : Note on the Importance of Imperfect-Gas Effects and 

Variation of Heat Capacities on the Isentropic Flow of Gases. NACA 
RM L8jl4, 1948. 

10. Baradell , Donald L. : Experimental Verification of Boundary-Layer Corrections 

in Hypersonic Nozzles. Jour. Aero/Space Sci. (Readers' Forum), vol. 26, 

no. 7; July 1959; PP- 454-455- 

11. McCracken, Daniel D. : A Guide to FORTRAN Programming. John Wiley & Sons, 

Inc., c. 1961 . 

12. Ferri, Antonio: Elements of Aerodynamics of Supersonic Flows. The MacMillian 

Co., 19^9- 

13. Beckwith, Ivan E., Ridyard, Herbert W. , and Cromer, Nancy: The Aerodynamic 

Design of High Mach Number Nozzles Utilizing Axisymmetric Flow With Appli- 
cation to a Nozzle of Square Test Section. NACA TN 2711, 1952- 


51 


l4. Persh, Jerome: A Theoretical Investigation of Turbulent Boundary Layer Flow 

With Heat Transfer at Supersonic and Hypersonic Speeds. NAVORD Rep. 3854, 
U.S. Naval Ord. Lab. (White Oak, Md.), May 19, 1955 . 

15* Cohen, Nathaniel B. : A Method for Computing Turbulent Heat Transfer in the 

Presence of a Streamwise Pressure Gradient for Bodies in High-Speed Flow. 
NASA MEMO 1-2-59L, 1959- 

l6. Bartz, D. R. : A Simple Equation for Rapid Estimation of Rocket Nozzle Con- 
vective Heat Transfer Coefficients. Jet Propulsion (Tech. Notes), vol. 27, 
no. 1, Jan. 1957> PP* 49-51. 


52 


NASA- Langley, 1963 L~3l80 


